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ABSTRACT 



This thesis describes the modification of six computer programs on the 
Micro VAX/2000/CAD/CAE workstation. Three of the programs, NEW.DOUBLE, 
NEW.PANEL, and NEW_VOR, were originally transferred to the Aeronautical 
Engineering VAX System Server by LCDR John Campbell. Two of the 
programs (SUB and SUPER), both vortex lattice method programs, were placed 
in the VAX system by Mr. Rich Margason of the Langley Research Center. 
None of the above five programs had any graphics facility. The sixth program, 
a viscous interaction program was transferred/adapted to the VAX system by 
the author of this report. Extensive modifications were subsequently made to 
these programs to enhance their user interface. In addition, all the programs 
have been adapted to provide interactive graphical/printed output. Further- 
more, program NEW_DOUBLE was modified to accept any arbitrary symmetri- 
cal shaped body. Lastly, NEW.PANEL was altered to interface with a viscous 
interaction effects program in which the boundary layer characteristics are 
determined. All user inputs in NEW_DOUBLE, NEW.PANEL and NEW.VOR were 
backed up with interactive checking routines. The programs are intended to 
be used by aeronautics/astronautics engineering students in basic and ad- 
vanced courses in aerodynamics. 
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THESIS DISCLAMER 

The reader is cautioned that computer programs developed in this 
research may not have been exercised for all cases of interest. While every 
effort has been made, within the time available, to ensure that the programs 
are free of computational and logic errors, they cannot be considered vali- 
dated. Any application of these programs without additional verification is at 
the risk of the user. 
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I. INTRODUCTION 



Incorporated in this thesis are six FORTRAN programs which have been 
extensively modified to improve their user interface and to enhance their out- 
put capabilities. Three of the programs, NEW_PANEL, NEW_DOUBLE, and 
NEW_VOR were originally transferred to the Aeronautical Engineering 
CAD/CAE Lab by LCDR John Campbell. Two of the programs, SUB and SUPER, 
vortex lattice programs, were developed by NASA AMES and transferred to the 
CAD/CAE Lab by Mr. Rich Margason. The sixth program, written by Dr. 
Cebeci, was transferred/adapted by the author of this thesis. 

The main focus of this thesis was to adapt the above mentioned pro- 
grams with a graphic facility. The graphics base program was developed by 
Mr. Dave Marco of the Mechanical Engineering Department, Naval Post- 
graduate School. The base program is simply a compilation of FORTRAN sub- 
routines (similar to the popular graphics program DISSPLA) which can be 
called to effect 2-D graphics. As such, the graphics produced by the aforemen- 
tioned Computational Fluid Dynamics (CFD) programs are limited to two 
dimensions. All of the plots generated by the respective program adaptations 
are effected interactively with little or no user input. 

Additionally, the respective program modifications were primarily in- 
tended to enhance the user interface with the CFD programs and, in a sense, 
effectively streamline their use. Specifically, adaptations such as input check- 
ing, interactive menus, automatic sorting routines, and backup output data 
file creation were integrated into each program. To supplement this thesis 
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objective, a concise User's Manual was created to provide the aeronautical 
engineering student a program reference guide. The user's manual, in fact, 
details all functional aspects of these programs and provides computational 
examples supplemented with graphical/tabular results. The User's Manual 
was geared to the student with little or no experience with the VAX system. 

Furthermore, the programs NEW_DOUBLE and NEW_PANEL each received 
specific adaptations which were not originally within the scope of this thesis. 
Specifically, the NEW_DOUBLE program, a line doublet distribution program, 
was adapted to input any symmetrical body. Originally, the NEW_DOUBLE 
program could only consider an elliptic or a one-family symmetrical airfoil- 
like shape; this family is described by the equation; 

Y(x) = A (c-x) 

The NEW_PANEL program was adapted to interactively process a coefficient of 
pressure (Cp) distribution to determine boimdary layer characteristics. The 
source of the Cp distribution to be analyzed can be either from the 
NEW_PANEL program itself or an arbitrary Cp distribution which can be 
entered from an input data file. 

Thesis results and recommendations for future work are given. 



2 



II. PROGRAM ADAPTATION 



A. INTRODUCTION 

Computer programs can always be improved or enhanced. New pro- 
gramming techniques coupled with improved programming languages provide 
a limitless number of programming innovations which can be initiated. Addi- 
tionally, as technological advances occur in the realm of computer hardware, 
it is without question that software development will also be enhanced. 

The first major modification made to all the programs considered in this 
thesis was to totally restructure and streamline each program to facilitate 
editing and compiling. As each program was modified, it became excessively 
large; the main program coupled with all of its subroutines to include the new 
graphics subroutines invariably exceeded the buffer size of the VAX 2000 
Workstation. The base program and each of its subroutines were placed in a 
separate FORTRAN file. Each FORTRAN subroutine was then compiled as a 
single entity to create its own object file. The object files of the subroutines 
were then consolidated into a library file specifically named to support its 
base program. Prior to running a particular program, the base program object 
file, the library of subroutine object files £md the object file for the graphics 
program were linked together to create a single executable file. 

All of the programs, with the exception of the viscous effects program, 
were adapted to output 2-D plots. As mentioned in the introduction, the 
graphics program was created by Mr. Dave Marco and, in its present state, is 
called DLIB. Again, DLIB is a compilation of graphics subroutines. The version 
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of DLIB used in this thesis is called DISL and represents a small alteration of 
the original to enhance its interface with the graphics subroutines. Specific- 
ally, the DLIB requirement to enter the command of "CONTINUE" once a 
graphics image is presented on the screen was deleted. The original intent of 
this requirement was to ensure that the graphics image would not be erased 
from the screen vmtil the user decided to continue. However, the graphics 
image created by a FORTRAN subroutine will not disappear imtil the program 
execution is terminated. Thus, this command requirement (CONTINUE) was 
not needed and presented a point of confusion for program users. 

The procedure to effect the graphical plots in each program differed to a 
degree. The principal difference lay in the method in which data arrays were 
read into the graphics subroutine. However, the use of backup data files was 
common to all programs. These files not only facilitated the development of 
the graphics but also acted as a data checking vehicle. The use of common 
blocks to transfer data arrays between subroutines was kept to a minimum 
for simplicity’s sake. Automatic scaling of data is not standardized across all 
of the CFD programs. Varying techniques were needed to preclude data array 
distortion. However, each graphics subroutine contains a FORTRAN "CALL" 
statement to effect the automatic scaling routine which determined the 
maximum and minimum values of the data arrays to be plotted. 

B. PROGRAM NEW.DOUBLE ADAPTATIONS 

The purpose of the NEW_DOUBLE program is to determine the piecewise 
constant doublet strength, m(t), for a line doublet distribution of an elliptic or 
s)mametrical airfoil-like shapes at zero angle of attack. The points tj, repre- 
sent the location of the doublets along the chord or line of symmetry. They are 
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concentrated near the ends of the distribution, using a cosine spacing method, 
where the variation of the doublet strength is expected to be most rapid. The 
point tj corresponds to Xg and tjvj corresponds to the endpoint xf. 

The stream function can be calculated from the doublet strength distri- 
bution. From the stream fimction, the velocity components and the pressure 
coefficients are then calculated. The surface shape is defined by y= Y(x) and 
the solution must satisfy zero velocity conditions at the leading and trailing 
edge stagnation points. 

In addition to adding graphics subroutines to this program, 
NEW_DOUBLE was adapted to analyze any symmetric shape. The user is first 
required to interactively enter the respective data points for the top portion of 
the s 3 Tnmetric shape. Once all of the points have been entered, the program 
will allow the user to correct any mistakes he or she may have made while 
entering the data. Using a spline routine, the intermediate points along the 
symmetric shape can be obtained readily to facilitate program processing 
[Ref. 1]. In brief, the spline routine created a continuous function between 
each adjacent data point. 

The NEW.DOUBLE graphics subroutines (GRAPHl, GRAPH2, GRAPHS) 
presented several unique features. First of all, the automatic scaling routines 
are specific to the particular type of plot to be created. Specifically, three 
automatic scaling routines were created: FIX, SCALER, SCALER2. Each of these 
routines determined the maximum and the minimum value of the specific 
array to be plotted. Another unique feature common to the NEW.DOUBLE 
graphics subroutines is that they were designed to read the data arrays to be 
plotted from dummy data files which were established in the computational 
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subroutines. This technique facilitated the data checking capability of the 
program and ensured accurate plots. Lastly, all of the plots were created to 
produce explanatory remarks to enhance the user's capability to relate the 
program inputs to the graphical outputs. The graphics subroutines included 
common blocks which enabled the plots to display user interactive inputs 
such as the thickness ratio, the maximum thickness, and the number of inter- 
vals specified. Typical plots obtained from the NEW_DOUBLE program are pro- 
vided in Appendix A. [Ref. 2 and 3] 

C. PROGRAM NEW.PANEL ADAPTATIONS 

The purpose of the original PANEL program was to provide an analysis of 
the aerodynamics of NACA four-digit airfoils and airfoils of the NACA 230XX 
family using the source panel method. The program has been modified to ac- 
cept arbitrary airfoil surface coordinate input and is limited to single-element 
airfoils. The solution is determined for conditions of incompressible and invis- 
cid uniform free-stream flow. The very small coefficient of drag provided in 
the results is due to numerical roimd-off error. Furthermore, NEW_PANEL has 
also been adapted to analyze viscous effects. When considering the viscous 
analysis loop of the program, it is important to understand that the Cebeci 
program adaptation is sensitive to flow separation on the airfoil. Boundary 
layer thickness and other botmdary layer characteristics are computed and 
outputted into a tabular format. 

The most dramatic modification made to the NEW_PANEL program was 
the adaptation of the program to consider viscous effects. The first step in 
making this modification was to transfer the Cebeci program to a CAD/CAE lab 
account. The original version of this program was provided by Dr. M. F. 
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Platzer, Aeronautics/Astronautics Department, Naval Postgraduate School 
(NPS). The program, written in FORTRAN, had already been adapted for an 
IBM PC but the user interface with this program was extremely poor and for- 
mal instructions on its use did not exist. The program was subsequently 
memually transferred to a VAX lab account. After an inordinate amount of 
error checking, the program was validated against a report offered by Dr. 
Platzer. The Cebeci program was then modified to enhance its user interface 
by incorporating interactive input requests to include printing options and 
input source selection. In addition, a common “bubble-sort” FORTRAN routine 
was added to the Cebeci program to automatically determine the stagnation 
point on the airfoil [Ref. 4]. The original version of this program required the 
user to specify this point. Furthermore, the user is required to specify the 
point at which laminar-t\irbulent transition occurs on both the top and the 
bottom of the airfoil as well as the stream-flow Reynolds number. Finally, the 
Cebeci program was fully integrated with the NEW_PANEL program. The Cp 
distribution created by the NEW_PANEL method was then interactively sorted, 
scaled, and inputted into the Cebeci program. This program, as noted above, 
then computed and outputted the respective boundary-layer characteristics. 
In addition to the Cp distribution created by NEW_PANEL program, the user 
can also enter any arbitrary Cp distribution from a data file called BL2D.DAT. 
This last option allows the user to in effect, conduct viscous effects calcula- 
tions while not being limited by the program restrictions of NEW_PANEL. 

The NEW_PANEL graphics subroutines (GRAFl, GRAF2) presented several 
unique features. First of all, the automatic scaling routines are specific to the 
particular type of plot to be created. Specifically, two automatic scaling 
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routines were created: FORMl, FORM2. Each of these routines determined the 
maximum and the minimum value of the specific array to be plotted. Like 
NEW_DOUBLE, NEW_PANEL graphics subroutines were designed to read the 
data arrays to be plotted from dummy data files which were established in 
the computational subroutines. Again, all of the plots were created to produce 
explanatory remarks to enhance the user’s capability to relate program in- 
puts to graphical outputs. The graphics subroutines included common blocks 
which enabled the plots to display user interactive inputs such as the number 
of panels, the angle of attack, and the NACA airfoil number. Another unique 
feature of the NEW_PANEL graphics subroutines is the addition of the capabil- 
ity to produce two graphs which were not within the scope of the original 
NEW_PANEL program. Specifically, the relationships of Cm c/4 versus angle of 
attack and Cl versus angle of attack can be plotted. This adaptation was real- 
ized by causing the NEW_PANEL program to perform the NEW_PANEL analysis 
at 2 degree increments in angle of attack from -8 degrees to 16 degrees. Typi- 
cal plots obtained from the NEW_PANEL program are provided in Appendix B. 
[Refs. 2 and 3] 

D. PROGRAM NEW.VOR ADAPTATIONS 

The purpose of the NEW_VOR progreim is to provide an application of the 
vortex lattice method for the determination of the lift distribution of a flat 
rectangular wing. This method is based on a distribution of discrete horseshoe 
vortices over a wing surface that has been divided into a finite number of 
panels. A system of linear equations is developed for the vortex strengths on 
the panels and solved by matrix methods. 



8 



In addition to adapting the NEW_VOR program for printing/graphics 
options, this program was also modified with a unique subroutine to effect the 
automatic scaling ftmction. Rather than creating a separate/unique scaling 
subroutine for each graphical output, a single subroutine (called MAXMIN) 
was developed to sort the designated array. The MAXMIN subroutine was de- 
signed to output the maximum and minimum value in the array, and the 
particular array in ascending order. In that the array to be plotted was out- 
putted in ascending order, it was necessary to establish a dummy array in 
each respective graphics subroutine which would be plotted. Otherwise, the 
plots would invariably ascend from left to right. Again, each plot was adapted 
to present user inputs. Specifically, the values for aspect ratio and angle of 
attack are displayed. Typical plots obtained from the NEW_VOR program are 
provided in Appendix C. [Ref. 2] 

E. PROGRAM "SUB" ADAPTATIONS 

The SUB program has been adapted from a National Aeronautics and 
Space Administration (NASA) FORTRAN program which has been used consid- 
erably at the Langley Research Center and in industry. The results have 
shown good correlation with experimental results. SUB has subsequently been 
revised to enhance it’s ease of use and its ability to present accurate graphical 
results. This particular program has also imdergone extensive student evalu- 
ations. An AE 2035 class of 14 students thoroughly tested and evaluated the 
majority of the functions which this particular program offers. As a result of 
their findings, humerous modifications were made to the program SUB as will 
be detailed below. 
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The purpose of the SUB program is to estimate the subsonic aerodynamic 
characteristics of complex planforms. The program represents a lifting plan- 
form with a vortex lattice. A relatively complex plemform may be analyzed 
using up to 24 hne segments on a semispan. Additionally, these line segments 
may have an outboard variable-sweep panel or they may have several dihed- 
ral 2 mgles across the span. Furthermore, two planforms may be used together 
to represent a combination of wings and tails or wing, bodies, and tails. 

The SUB graphics subroutines (GRAPHl, GRAPH2, GRAPHS) present several 
unique features. First of all, automatic scaling is again effected by the 
MAXMIN routine. Secondly, dummy data files were established in the compu- 
tational subroutines and subsequently read in each graphics subroutine. The 
use of common blocks was kept to a minimum. The coefficient of pressure 
data provided by the SUB program lends itself readily to 3-D graphics. How- 
ever, in the absence of a 3-D graphics program in the CAD/CAE Lab, the pro- 
gram was modified to locate the data at the user specified planform position. 
Specifically, a sorting routine was developed to allow the user to specify a 
particular spanwise position on the planform to analyze the Cp distribution 
across the chord of the planform. In order to create this sorting routine, it was 
necessary to adapt the data output to the finite difference nodal network. 
This was simply done by recdizing the constant spacing distances between the 
nodal points (stations). T 3 q>ical plots obtained from the SUB program are pro- 
vided in Appendix D. [Ref 5] 

Lastly, the SUB program was modified to provide the user the opportu- 
nity to copy the output data file into an alternate data file so that his or her 
results would be saved for further gmalysis. Subsequent runs of the program 
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could then be made without losing the results already determined. This modi- 
fication was effected through interactively allowing the user to select an 
alternate data file from a list of four files. 

F. PROGRAM •’SUPER" ADAPTATIONS 

The SUPER program has also been adapted from a National Aeronautics 
and Space Administration (NASA) FORTRAN program and has been used con- 
siderably at the Langley Research Center. The use of this program is confined 
to the supersonic flow regime. In addition, the linearized supersonic lifting 
surface theory, used in this program, applies to wings having negligible thick- 
ness . SUPER has subsequently been revised to enhance it’s ease of use and its 
ability to present accurate graphical plots. These graphical representations 
have been verified with NASA reports as referenced below. 

The purpose of the SUPER program is to estimate the supersonic aero- 
dynamic characteristics of complex planforms. Linearized supersonic lifting 
surface theory is employed to calculate the aerodynamic characteristics of a 
warped wing of arbitrary planform. The program calculates lifting pressure 
distribution for the warped wing at fixed attitude and the pressure distribu- 
tion (per degree angle of attack) for a corresponding flat wing. These two 
pressure distributions are combined by superposition principles and inte- 
grated over the wing surface to obtain the variation of aerodynamic charac- 
teristics with changes in angle of attack. 

Similar to the case of program SUB, complex sorting routines were devel- 
oped to allow the user to specify the respective chordwise or spanwise position 
on the supersonic planform which would be analyzed for plotting purposes. 
The coefficient of pressure data provided by the SUPER program also lends 
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itself quite readily to 3-D graphics. Again, in the absence of a 3-D graphics 
program, the program was modified to locate the data at the user specified 
planform position. In order to create these sorting routines it was necessary to 
adapt the data output to the finite difference nodal network. Like program 
SUB, this was done by realizing the constant spacing distances between the 
nodal points and subsequently sorting the data accordingly to isolate the re- 
quested data. Typical plots obtained from the SUPER program are provided in 
Appendix E. [Refs. 6 and 7] 

Another modification made to the SUPER program concerns its output 
data file. The output data file for the SUPER program is extremely long. The 
great length of this file was of negligible utility to the common user. A new 
output data file was created within the text of the program which simply 
outputted the input data and the aerodynamic results. The tabular coefficient 
of pressure data was not incorporated into the output file. However, the 
full output data file with complete Cp data is still written to a file called 
’’OUTER.DAT". The abbreviated output file (OUTFILE.DAT) greatly facilitates 
the printing of the output file for the user. 

Lastly, like SUB, the SUPER program was modified to provide the user 
the opportunity to copy the output data file into an alternate data file so that 
his or her results would be saved for further analysis. The user is given the 
opportunity to interactively select an alternate data file name in which he or 
she can store their computational results. 
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III. SOLUTION FOR THE TWO-DIMENSIONAL 

INCOMPRESSIBLE LAMINAR AND TURBULENT 
BOUNDARY LAYER PROBLEM 

A. INTRODUCTION 

This section relates the numerical methods employed to solve the two 
dimensional incompressible laminar and turbulent boundary layer problem 
as conceived by Dr. T. Cebeci and Dr. H. B. Keller. As discussed earlier, this 
particular boundary-layer solution method was modified and imbedded into 
the NEW_PANEL program. The intent of this section is to provide a brief syn- 
opsis of their problem solution, not a detailed account. The development of 
the specific theoretical basis/computer code development of the Cebeci pro- 
gram is not within the scope of this thesis [Ref. 8]. In order to use Dr. Cebeci's 
method, it is necessary to input the potential flow solution over a section 
shape. In particular, the Cp distribution or the velocity distribution is 
required. Such information is obtained quite readily through the execution of 
the NEW_PANEL program. In fact, the Cp distribution is interactively sorted 
and inputted to the Cebeci program upon the user's request. In addition, one 
of the functional capabilities of the NEW_PANEL program is to input an 
arbitrary velocity distribution. Furthermore, the Cebeci program version 
provided by Dr. M. F. Platzer was further revised to determine the coeffi- 
cients for skin friction drag and form drag from the computed boundary-layer 
characteristics. This additional capability was transferred from the original 
version which is currently available for use on the IBM mainframe at the 
Naval Postgraduate School, account 4632P. 
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B. NUMERICAL SOLUTION BASIS 



This program uses a finite-difference method to solve the partial differ- 
ential equation obtained by using the Falkner-Skan transformation of the 
general boundary layer equations. Both laminar and turbulent flows may be 
analyzed in that an eddy-viscosity concept has been incorporated into the 
program which allows the momentum equation for turbulent flows to be writ- 
ten in the same form as a laminar flow. Dr. Cebeci's method is valid except 
upon the evolution of flow separation. The boundary layer separation point 
corresponds to the vanishing of the wall shear force at that point. Dr. Cebeci 
[Ref. 8] states, "if the wall shear vanishes at some x-location during the 
solution procedure, the solutions break down and convergence cannot be 
obtained. This is sometimes referred to as the singular behavior of the wall 
shear close to the separation point." Close inspection of the boundary layer 
results provided by the NEW_PANEL program is advised in order to ensure 
that the results are in fact valid. Extremely large values of displacement or 
momentum thickness indicate flow separation on the shape being analyzed. 
As an additional note, the program is limited to two dimensions in that 
negligible transverse curvature has been assumed. 

C. COMPUTER PROGRAM NUMERICAL SOLUTION 

There exist several methods to solve the boundary-layer equations. The 
finite difference method used in this program was developed by Dr. H. B. 
Keller [Ref. 9]. Keller's box method has been used extensively to solve the 
boundary-layer equations. The first requirement to be effected before the 
Keller method can be employed is to rewrite the governing equations as a 
first order system. The resulting first-order equations are subsequently 
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approximated on an arbitrary rectangular net. The finite-difierence equations 
evolve from "centered-difference" derivatives and averaged at the midpoints 
of the net rectangle. Figure 1 represents the orientation of the net rectangle. 
The respective nodal points are determined by: 

Eq— E jj = n = 1,2,3...N 

nQ=0, nj = nj.i +hj, j - 1,2, 3... J 




Figure 1. Rectangular Net Orientation, Keller's Box 



Solutions of the finite difference equations yield a truncation error of the 
second order. The difference equations are subsequently linearized by New- 
ton's method [Ref. 9]. Finally, the equations are solved by a block-elimination 
method [Ref. 8]. 

The computer program has been broken down into several separate sub- 
routine programs labeled as CIB, COEF, BL, EDDY, SOLV3, OUTPUT, and DRAG. 
Subroutine CIB is called from the NEW_PANEL main program which in turn 
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calls the remaining subroutines in order to determine the requisite boundary- 
layer characteristics. The flow transition point Gaminar to turbulent) is inter- 
actively inputted by the user for both the upper and lower surfaces of the 
airfoil. The user is also prompted to enter the chord-based Reynolds number. 
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IV. USER'S MANUAL 



In order to facilitate the use of the programs addressed in this thesis, a 
User's Manual was created which expanded upon the one made by LCDR 
John Campbell [Ref. 1]. Appendices A through E of this thesis contain the 
text portions of the User’s Manual as well as representative graphical out- 
puts. The User's Manual in its final form is approximately 150 pages long, 
excluding section and sample problem dividers. The bulk of the User's 
Manual consists of the output data files, the input data files (if appropriate) 
and the graphical outputs for each sample problem referenced in the User's 
Manual. These output files/plots served as the primary basis for validating 
the graphical plots obtained by each respective program. As an additional 
note, the User's Manual was not included in this thesis in its entirety due to 
its extreme length. 
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V. PROGRAM COMPUTER CODES 



Appendices F through J contain the complete source codes of the pro- 
grsuns discussed in this thesis. These codes are in their final form. However, 
it should be noted that each subroutine/main program has been written and 
presented as a stand-alone FORTRAN program in that each program can be 
compiled individually. As noted earlier, once each subroutine was compiled, 
an "object" file was created by the computer. This "object" file was then placed 
in its respective program library (DOUBLIB, PANLIB, VORLIB, SUBLIB, SUPLIB). 
Prior to running the particiilar program, the library of "object" files was 
linked with the object file of the main program (NEW_DOUBLE, NEW_ PANEL, 
NEW_VOR, SUB, SUPER) and with the object file of the graphics program (DISL). 
Linking the files in this manner created a single executable file for each 
program. 

There exists several lines of FORTRAN code in each individual program 
which are currently not executable (comment lines). Some of these comment 
lines will facilitate future modifications; others represent routines which were 
specifically incorporated for data checking; and the remainder are simply 
comments to clarify the executable statements. Rather than deleting these 
lines as being extraneous, they were left in the program to facilitate future 
modifications/program maintenance. These routines are marked appropri- 
ately within the program. 
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VI. RESULTS AND RECOMMENDATIONS 



The objectives of this thesis, as originally conceived, have been realized. 
Five FORTRAN programs (NEW.DOUBLE, NEW.PANEL, NEW.VOR, SUB, SUPER) 
have been modified to interactively supply graphical representation of the re- 
spective computational results. In fact, SUB and SUPER were added to the list 
of programs to be modified well into the thesis research process. In addition, 
NEW_DOUBLE can now analyze any symmetrical shape. Data checking rou- 
tines were also added to NEW_DOUBLE to enhance data input procedures. 
Furthermore, the user interface capabilities of each program were signifi- 
cantly improved, especially for SUB and SUPER. Lastly, each program was 
adapted to provide the user the capability to interactively print the computa- 
tional results or the plots developed. 

To enhance the utility of these programs, a concise User’s Manual was 
developed. This manual fully describes each program to include program and 
input restrictions. In addition, numerous sample problems were integrated 
into the manual to demonstrate to the user the various capabilities of each 
program. For each sample problem, detailed instructions are given on how to 
use the program properly. Furthermore, the output data files and graphical 
plots created by each sample problem are also included in the User's Manual. 

The validation of the graphical results was achieved through several 
sources; in particular, the books by Ira H. Abbott and A.E. VonDoenhoff 
[Ref. 10] and by John D. Anderson, Jr. [Ref. 11], were used to check the plots 
generated by NEW.DOUBLE and NEW.PANEL. LCDR J.A. Campbell's thesis 
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results were used to validate the NEW_VOR plots. In order to check the SUB 
and SUPER graphs, the NASA publications detailing each respective program 
were used [Refs. 5, 6, and 7]. In all cases, the graphical representations 
produced by these programs were qualitatively and quantitatively correct. No 
arbitrary ac^justments were made to the graphics subroutines to "fit" the data 
to the respective validating source. 

Modifications and further adjustments can always be made to a com- 
puter program to either enhance or expand its capabilities. As stated earlier, 
the tabvJar output of data in SUB and SUPER readily lends itself to 3-D graph- 
ics display. At such time that the Aeronautics/Astronautics Department 
CAD/CAE lab receives a 3-D general graphics package, such as DISSPLA, SUB 
and SUPER can easily be adapted to produce 3-D plots. The graphics subrou- 
tines, as they are currently written, use call statements identical to those 
used with a DISSPLA package. Furthermore, the data generation required for 
the respective 3-D plots has already been programmed into the graphics sub- 
routines. An additional modification would be to adapt the Cebeci program 
output to produce graphical results. Furthermore, the Cebeci program could 
also be modified to solve a variety of problems including 2-D flows with heat 
and mass transfer, slot injection as well as axisymmetric flows. Lastly, the 
programs SUB and SUPER could be adapted to interactively accept the data 
inputs from the console rather than requiring the user to create an input data 
file. However, the inputs to both programs can be rather long and detailed in 
the analysis of a complex planform. 
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I. INTRODUCTION 



The purpose of the NEW_DOUBLE program is to determine the piece- 
wise constant doublet strength m(t) for a line doublet distribution of an ellip- 
tic or airfoil-like shape at zero angle of attack. The points tj, represent the 

location of the doublets along the chord or line of symmetry. They are concen- 
trated near the ends of the distribution, using a cosine spacing method, where 
the variation of the doublet strength is expected to be most rapid. The point tj 
corresponds to Xg and tjsf corresponds to the endpoint xf. 

The stream fimction can be calculated from the doublet strength distri- 
bution. From the stream fimction, the velocity components and the pressure 
coefficients may be calculated. The surface shape is defined by y=Y(x) and the 
solution must satisfy the zero velocity conditions at the leading and trailing 
edge stagnation points. 

II. ASSUMPTIONS AND LIMITATIONS 

The approach taken to develop this method of solution assumes that the 
doublet strength functions are both piecewise-constant along the chord. It is 
also important to remember that this solution is valid for incompressible and 
inviscid uniform freestream flow. Since the bodies under investigation are 
(two dimensional) symmetrical and at zero angle of attack, there is no lift nor 
induced drag produced. In addition, there is no drag since we are considering 
an inviscid fluid and no separation is allowed for. 
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III. INPUT DESCRIPTION 



There are very few input values required for this simple program. Their 
description and program variable names are listed below. 



NTYPE — ^Type of body shape; elliptic, a single-family airfoil-like, given by 



Y(x) = A 




(c-x) 



or symmetric. 



TAU — Thickness ratio. (Maximum thickness/chord) 

XMAXY — Chordwise location of the point of maximum thickness. (Airfoil only) 
N — Number of intervals. 2 < N < 100 
XS — Doublet distribution starting point. 

XF — Doublet distribution ending point. 

NXTOL — Exponent value used to generate the convergence criterion XTOL. 
NFTOL — Exponent value used to generate the convergence criterion FTOL. 
XTOL — X location tolerance. 

FTOL — Function tolerance. 

X-X Coordinate of the symmetric shape airfoil surface. 

Y-Y Coordinate of the symmetric shape airfoil surface. 



IV. SAMPLE PROBLEMS 

A few sample problems will illustrate the use of the NEW_DOUBLE 
program. The first problem will use an ellipse of thickness ratio 0.1. The 
second problem will analyze an airfoil-like shape with a thickness ratio of 
0.12 and a chordwise location of maximum thickness of 0.30. Finally, the 
third problem will analyze a symmetric shape. 
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V. STARTING THE PROGRAM 



Begin with the screen showing the DCL prompt, which looks like this: 

$ 



Next, ensure that the program is in your directory by typing: 

DIR [Return] 

and viewing the files for NEW_DOUBLE.EXE. 

To run the program, type: 

RUN NEW_DOUBLE [Return] 

The program will start and the screen should look similar to that shown 
in Figure 2. 



PROGRAM NEW_DOUBLE; VERSION 3 : 4 OCTOBER 89 



DOUBLET DISTRIBUTION METHOD IS USED TO DETERMINE 
INCOMPRESSIBLE FLOW AROUND AN ELLIPSE . SYMMETRICAL 
AIRFOIL OR ARBITRARY-SYMMETRIC SHAPE AT ZERO ANGLE OF ATTACK 

PROGRAM ASSUMES A NONDIMENSIONAL CHORD, THAT IS, THE 
VALID RANGE OF X IS FROM 0 TO 1. 

ENTER TYPE OF BODY SHAPE DESIRED: 

1) ELLIPTIC 

2) SYMMETRICAL AIRFOIL-LIKE OR 

3) ARBITRARY SYMMETRIC SHAPE 

ENTER 1,2, OR 3. 

NOTE THAT OPTION 3 WILL REQUIRE MANUALLY INPUTTING DATA 
POINTS FOR THE UPPER SIDE OF THE RESPECTIVE BODY 



Figure 2. Initial Screen for Program NEW_DOUBLE 
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VI. SAMPLE GRAPHICAL OUTPUTS 
A. SAMPLE PROBLEM ONE 

For the elliptic case, respond to the initial screen request by 

entering: 

1 (Return) 

Respond to the request for the thickness ratio by entering: 

0.1 (Return) 

Now enter the number of intervals you desire the doublet distri- 
bution to have by entering: 

10 [Return] 

The screen should now look like that shown in Figure 3. 



WHICH METHOD DO YOU WISH TO USE TO DETERMINE THE DOUBLET 
DISTRIBUTION ENDPOINTS? 

1) PROGRAM INTERVAL-HALVING SUBROUTINE TO ITERATE 

2) MANUAL ITERATION BY THE USER 

3) RETURN TO START 

4) EXIT PROGRAM 
ENTER 1,2,3 OR 4 



Figure 3. Endpoint Determination Method Selection Screen 
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Respond to the question by entering: 

1 [Return] 

If you should desire to enter your own values, enter 2. 

The next values you will be reqmred to enter are for the X location 
tolerance and the stagnation point velocity function tolerance. It is recom- 
mended that values of lOE-6 (0.000001) be used. The maximum number of 
iterations should be set at a value of at least 20 when using such small 
tolerances. Additionally, if you desire to use, for example, 10 intervals, you 
should use lOE-4 so as to achieve a small velocity magnitude at the stagna- 
tion points. 

The output parameter entry has only to do with the interval halving 
subroutine. Unless you are having problems with the program or are inter- 
ested in the convergence of the solution, it is recommended that this value be 
set to zero (0). 

Following entry of the output parameter, the program begins the 
solution process. It returns with UO and Ul, the values for the X velocity 
component at the leading and trailing stagnation points respectively and the 
values for XS and XF, the beginning and ending points of the line doublet 
distribution. If the values for UO and Ul are sufficiently close to zero, say 
less than lOE-3 (0.001), then enter: 

Y [Return] 

If you desire more accuracy, enter: 

N [Return] 
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and then reenter the tolerance and maximum iteration values. Responding 
with a (Y) will cause the progrsun to proceed to the output stage. Values will 
be printed to the screen and to the following data files: 

DUBLET.DAT : DOUBLET STRENGTH DISTRIBUTION 

SHAPE.DAT : BODY SURFACE COORDINATES 
PRESSURE.DAT : SURFACE PRESSURE DISTRIBUTION 

You will be asked for the number of pressure coefficient output 
points you desire. This number is independent of the number of intervals of 
the line doublet distribution. It affects only the number of output data points 
and not the accuracy of the solution. After entering the number of Cp output 
points, pressure distribution data will be displayed to your screen. The 
program now asks if you want to print the results (Y/N). Enter your response 
and select the respective file which you want to print from a tabulated listing. 
However, be aware that you must have already logged onto the KELLY 
terminal to print anything, or be at a terminal which is connected to a 
printer. 

You will now be asked if you want to graph the results (Y/N). If you 
respond affirmatively, the screen will look similar to Figure 4. 

Once you have selected your plotting option and the respective plot 

has appeared on your screen (on the KELLY terminal screen if you are 

printing items) you will be asked if you would like a print of the plot (Y/N). 

Answer accordingly and continue with the program. 

» 

You will now be asked if you would like to make another run. Enter: 

1 [Return] 
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WHICH OF THE FOLLOWING DATA FILES 
DO YOU WANT TO GRAPH? 

1) DUBLET.DAT 

2) PRESSURE.DAT 

3) SHAPE.DAT 

4) NONE 

INPUT OPTION NO. (1 ,2,3 OR 4) 



Figure 4. Plotting Options Screen 



B. SAMPLE PROBLEM TWO 

Sample problem two will work through the airfoil-like shape case 
and the user will supply the values of XS and XF. The user may experiment 
with manual iteration, however to save space this sample will use previously 
determined satisfactory values of XS and XF for the initial guess. 

You should now be back at the initial screen and it should look like 
Figure 2. For the airfoil-like case enter: 

2 [Return] 

Respond to the request for the thickness ratio by entering: 

.12 [Return] 
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For the chordwise location of maximum thickness, enter: 

.30 [Return] 

Now enter the number of intervals you desire the doublet distribution to have 
by entering: 

10 [Return] 

The next step is to select the method for the determination of the 
endpoints for the doublet distribution. The screen should look like Figure 3. 
This time respond to the question by entering: 

2 [Return] 

For the doublet distribution starting point, XS, enter 

.0082129128 [Return] 

For the doublet distribution ending point, XF, enter 

.9994138 [Return] 

As with the previous example, the program now begins the solution 
process. It returns with UO and Ul, the values for the X velocity component 
at the stagnation points. It also echoes back the values entered for XS and 
XF. If the returned values for UO and Ul are sufficiently close to zero, then 
enter: 

Y [Return] 

This response will cause the program to proceed to the output stage. 
Values will be printed to the screen and to the data files. 
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Enter the number of pressure coefficient output points you desire. 
You are reminded that this number is independent of the number of intervals 
of the line doublet distribution and it does not affect the accuracy of the 
solution. 

Again, you will be afforded the opportunity to print and graph the 
results as in sample problems one. 

The program now asks if you want to make another run. Enter: 

1 [Return] 

C. SAMPLE PROBLEM THREE 

Sample Problem Three provides an example of arbitrary shape 
analysis. You should now be back at the initial screen and it should look like 
Figure 2. For the symmetric shape case enter: 

3 [Return] 

Once you have entered this response, your screen should look 
similar to Figure 5. 



HOW MANY UPPER PROFILE DATA POINTS DO 
YOU DESIRE? (ENTER A NUMBER BETWEEN 3 AND 100) 

BE AWARE THAT THE LEADING EDGE OF YOUR DESIRED 
SHAPE HAS BEEN PROGRAMMED TO BE AT THE ORIGIN 
AND THAT YOUR TRAILING EDGE IS AT (1,0). SCALE 
YOUR SHAPE/OBJECT ACCORDINGLY. 



Figure 5. Symmetric Shape Data Point Input Screen 
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Enter the number of points you wish to use to describe your 
symmetric shape. You will then be given the opportunity to enter each point. 
Once you have entered all of your surface points, the program will ask if you 
want to check your input data. You may then make any corrections as 
necessary. When you have finished correcting your data, enter N to the 
question asking you if you have any ‘input data corrections.’ The program will 
then proceed as described in example problems one and two. 

This completes the sample problems for the NEW_DOUBLE 
program. Representative graphical outputs created by these sample nms are 
listed in Figures 6 through 8. Since the bodies analyzed by this program are 
symmetrical with respect to the x axis, only the upper surface body shape 
coordinates and pressure coefficients are output. For this reason, the 
piecewise constant doublet strength M(I) is divided by two to indicate the 
portion affecting the upper surface. 
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X 

Figure 6. Doublet Strength Distribution 
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-0.312 -0.069 0.173 0.415 0.658 0.900 




X 



Figure 7. Cp Distribution 
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p. 000 0. 092 0. 184 0. 276 0. 368 0. 460 




X 

Figure 8. Airfoil Shape 
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I. INTRODUCTION 



The purpose of the NEW_PANEL program is to provide an analysis of 
the aerodynamics of NACA four-digit airfoils and airfoils of the NACA 230XX 
family using the panel method. This program has been modified to accept 
arbitrary airfoil surface coordinate input. NEW_PANEL has also been 
adapted to analyze viscous effects. 

II. ASSUMPTIONS AND LIMITATIONS 

This program is limited to single-element airfoils. The solution is deter- 
mined for conditions of incompressible and inviscid irrotational flow. The co- 
efficient of drag provided in the results is due to numerical round-off error. 
When considering the viscous analysis loop of the program, it is important 
that you understand that the Cebeci eddy-viscosity program adaptation is 
sensitive to flow separation on the airfoil. Boundary layer thickness and other 
boundary layer characteristics will be computed. It is advised that viscous 
analysis be limited to small angles of attack and to relatively slender airfoils. 

III. INPUT DESCRIPTION 

As with the NEW_DOUBLE program, there are very few input values 

required for this simple program. Their description and program variable 

names are listed below. 

NUPPER - Number of nodes on the upper surface. 

NLOWER - Number of nodes on the lower surface. 

X(1),Y(1) - Surface coordinates. 

These may be entered from the keyboard, from a data file, or from data 
statements. The program is capable of generating an approximation for air- 
foils of the NACA XXXX and 230XX series. 
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ALPHA — ^Angle of attack. (Angle between the chord and the freestream velocity.) 
RL — Chord Reynold’s Number 

XCTRI(l) - Flow transition point from laminar to turbulent flow on the top 
of the airfoil 

XCTRK2) - Flow transition point from laminar to turbulent flow on the 
bottom of the airfoil. 



IV. INPUT RESTRICTIONS 

The program, as written, is limited to 100 total surface nodes. This may 
be modified by changing the size of the arrays; however, only a very complex 
surface should require that many values to accurately define the surface. If 
that is the case, a more sophisticated program should be considered for the 
investigation. As mentioned above, the computer generated approximations 
to airfoil shapes are limited to the NACA XXXX and 230 XX series. The pro- 
gram will accept values for ALPHA up to 90 degrees, but the user is cau- 
tioned that since separation can exist for angles of attack as low as 5°— 10°, 
results for values above about 10° may be suspect. 

V. OUTPUT VARIABLES FOR VISCOUS RESULTS 

XC - Airfoil Coordinate (Abscissa) 

S - Station Location 
VW - Wall Shear 
CF - Coefficient of Friction 
DLS - Displacement Thickness 
THT - Momentum Thickness 
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VI. SAMPLE PROBLEMS 

A few sample problems will illustrate the use of the NEW_PANEL pro- 
gram. The first nm will be done using an approximation to a NAG A 0012 air- 
foil which is generated by the program using the information associated with 
each digit in the NACA number. The second run will analyze a NASA 
LS(1)-0013 airfoil using a set of data statements containing the airfoil surface 
coordinates. These statements have been inserted into the proper location in 
the program already. The last sample problem will re-analyze the LS(1)-0013 
airfoil but now viscous effects will be included. 

VII. STARTING THE PROGRAM 

Begin with the screen showing the DCL prompt, which looks like this; 

$ 

Next, ensure that the program is in your directory by typing: 

DIR [Return] 

and viewing the files for NEW_PANEL.EXE. To run the program, type: 

RUN NEW_PANEL [Return] 

The program will start gmd the screen should look similar to what is 

shown in Figure 9. 
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PROGRAM NEW_PANEL 

SMITH-HESS (DOUGLAS) PANEL METHOD FOR A SINGLE-ELEMENT 
LIFTING AIRFOIL IN TWO-DIMENSIONAL INCOMPRESSIBLE FLOW 

DO YOU WISH TO: 

1) USE AIRFOIL SURFACE COORDINATE DATA VALUES. 

2) HAVE COMPUTER GENERATE AN APPROXIMATION FOR 
NACA XXXX OR 230XX AIRFOIL SECTION. 

3) QUIT THE PROGRAM 
ENTER 1.2. OR 3 



Figure 9. Initial Screen for Program NEW_PANEL 



VIII. SAMPLE GRAPHICAL OUTPUTS 
A. SAMPLE PROBLEM ONE 

For the first case we will have the computer generate an approxi- 
mation for the shape of a NACA 0012 airfoil, consisting of 20 surface panels, 
using an algorithm contained in subroutine NACA45. The angle of attack of 
the onset flow will be six degrees. To use the approximation method, enter: 

2 [Return] 

Respond to the request for the ntomber of surface data points by entering: 

20 [Return] 
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Confirm the number of surface data points you desire by entering: 

1 [Return] 

Although the program will allow a different number of upper and lower sur- 
face data points, it is recommended that you try and keep them equal. An 
imequal nmnber of nodes yields trailing-edge panels of imequal length, which 
lowers the accuracy of the approximation of the Kutta condition. Respond to 
this question by entering: 

1 [Return] 

The next question asks for the NACA number of the airfoil you are consider- 
ing. For this case we will look at the NACA 0012, so enter: 

0012 [Return] 

The screen should now look like similar to Figure 10. 

The program is now ready to perform its calculations. The final 
piece of information required is the angle of attack, ALPHA. For this case, re- 
spond to the question by entering: 

6 [Return] 

Following entry of the angle of attack, the program begins the solution pro- 
cess. Values scroll up the screen and are simultaneously being written to the 
data files. You should now see a screen similar to the one shown in Figure 11. 

Should you select to print the resvilts, you will be given the option 
to print both of the data files or just the one you want. Once you have finished 
printing the results, you will be asked if you want to graph the results. 
Respond affirmatively and the screen should then look similar to Figure 12. 
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ENTER NUMBER OF SURFACE DATA POINTS DESIRED 

20 

NUMBER OF SURFACE DATA POINTS TO BE GERATED = 20 
IS THIS VALUE CORRECT? (YES=1. NO=2) 

1 

ARE THE NUMBER OF UPPER AND LOWER SURFACE 
DATA POINTS (NODES) EQUAL? (YES=1 . (NO=2) 

1 

INPUT NACA NUMBER, ANY FOUR DIGIT OR 230XX SERIES 

0012 

INPUT ALPHA IN DEGREES 



Figure 10. Screen Showing Data for Computer Generated Airfoil 



PROGRAM NEW_PANEL RESULTS HAVE BEEN WRITTEN TO FILES: 
PBODY.DAT : BODY SURFACE COORDINATES 
PPRESS.DAT : SURFACE PRESSURE DISTRIBUTION 
WOULD YOU LIKE TO PRINT THE RESULTS (Y/N)? 



Figure 11. Printing Option Screen 
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WHICH OF THE FOLLOWING DATA OUTPUTS 
DO YOU WANT TO PLOT? 

1 ) PPRESS.DAT (CP DISTRIBUTION) 

2) PBODY.DAT (AIRFOIL SHAPE) 

3) CL VS. ANGLE OF ATTACK 

& CM C/4 VS. ANGLE OF ATTACK 

4) NONE 

INPUT OPTION NO. (1 ,2,3 OR 4) 



Figure 12. Graphical Selection Screen 



Once the selected plot is displayed on your screen (screen KELLY 
if you are printing) you will be given the option of printing the plot. Again, 
you must have already used the "set host kelly" command to print items. If 
you elect not to print the graphical output you screen will again look similar 
to Figure 12. Selecting option 4 (NONE) will exit you from the graphing loop. 
You will now be asked to analyze the viscous effects for the airfoil. Respond 
negatively by entering: 

N [Return] 

A new screen will be presented emd the program now asks if you want to 
make another run. Enter: 

1 [Return] 
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B. SAMPLE PROBLEM TWO 

This time the sample problem will examine a NASA LS (D-0013 
whose coordinates have been entered as data statements in the program. You 
should now be back at the initial screen and it should look like Figure 9. 
Since you will be using actual airfoil coordinate data values, enter; 

1 [Return] 

The screen shown in Figure 13 now presents you with the three 
choices available for entering the airfoil surface coordinate data values. You 
will be using the data statements, so enter: 

3 [Return] 



DO YOU WISH TO ENTER THE SURFACE COORDINATE VALUES; 

1) FROM A DATA FILE. 

2) FROM THE KEYBOARD. 

3) USING DATA STATEMENTS ALREADY ENTERED IN THE 
MAIN PROGRAM. “NOTE** THIS REQUIRES THAT 
PROGRAM BE MODIFIED IN ADVANCE BY MOVING DATA 
STATEMENTS TO THE CORRECT LOCATION. 

ENTER 1, 2. OR 3. (FOR PREVIOUS MENU ENTER 4) 



Figure 13. Menu for Surface Coordinate Data Entry Method 
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The number of data points has been entered via the data state- 
ments, therefore you are not asked that question for this case. For the angle 
of attack, again enter: 

6 [Return] 

As you saw in the previous example, values scroll up the screen. 
The program will again allow you to print or graph the respective results as 
before. Additionally, you will again be asked if you want to analyze viscous 
effects. Respond accordingly to exercise the required program options. 
Finally, the program will ask if you want to make another run. Enter: 

1 [Return] 

C. SAMPLE PROBLEM THREE 

As noted earlier, this sample problem will again analyze the 
LS(1)-0013 airfoil but with viscous effects. You should now be back at the ini- 
tigJ screen and it should look like Figure 9. Since you will be using airfoil sur- 
face coordinate data values enter: 

1 [Return] 

The screen should again look like Figure 13. Again enter the response 3 in 
that data statements will again be used. 

3 [Return] 

For the angle of attack response enter: 

0 [Return] 
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As you have seen in the two previous examples, values scroll up 
the screen. The program will again allow you to print and/or graph the 
respective results as before. When asked if you would like to analyze the 
viscous effects for this airfoil enter; 

Y [Return] 

The screen should now look similar to Figure 14. 

The first option (1 ) is used to input an arbitrary external velocity 
profile. The external velocity values at each respective point were obtained 
from the expression: SQRT(l-Cp). To input the Cp distribution just created 
for the LS(1)-0013 airfoil enter: 

2 [Return] 



VISCOUS BOUNDARY LAYER ANALYSIS 

*** INPUT DATA OPTION *** 

WHAT INPUT SOURCE WOULD YOU LIKE TO USE ? 

1) DATA FILE "BL2D.DAT" OR 

2) NEW_PANEL CP DISTRIBUTION JUST CREATED 

3) QUIT PROGRAM 

ENTER 1,2, OR 3 



Figure 14. Menu for Viscous Data Input Option 
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You will now be asked to enter the flow Reynold's Number. Enter; 

6000000 [Return] 

Now you will be asked to enter the respective nondimensionalized 
values of XCRIT(l) and XCRIT(2). Ageiin these values correspond to the point 
along the chord of the airfoil at which flow transition from laminar to turbu- 
lent occurs for the top and bottom of the airfoil, respectively. Enter .3 for both 
values. To avoid flow separation, these value should be greater than 0.15 for 
analysis at angles of attack in excess of approximately 5°. 

The program will now begin to process the input data and deter- 
mine the boundary layer characteristics. Upon completion of the computa- 
tions the screen should look similar to Figure 15. 

Remember that in order to print the results you must be logged 
onto a computer terminal which is connected to a printer. Enter the following 
command to print the boimdary layer results: 

Y [Return] 

Once you have finished the viscous flow analysis process, the program will 
again ask you if would like to make another run of NEW_PANEL. Enter the 
following command to exit the NEW_PANEL program: 

2 [Return] 

This completes the sample problems for the NEW_PANEL 
program. Representative graphical outputs created by these sample nms are 
shown in Figures 16 through 18. 
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READING THE DATA... 

INPUT OF DATA COMPLETE. 

BOUNDARY LAYER COMPUTATIONS IN PROGRESS... 
BOUNDARY LAYER COMPUTATIONS IN PROGRESS... 

THE BOUNDARY LAYER RESULTS HAVE BEEN 
WRITTEN TO FILE "BL2D.OUT" 

WOULD YOU LIKE TO PRINT THESE RESULTS ? 



Figure 15. Viscous Data Output File Option Screen 
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- 2.544 - 1.869 - 1.194 - 0.518 0.157 0.832 




X 



Figure 16 . Cp Distribution 
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-0.500 -0.300 -0.100 0.100 0.300 0.500 




X 

Figure 17. Body Shape 



49 



MOMENT <C/4) & LIFT COEFFICIENTS 
-1.461 -0.890 -0.319 0.252 0.823 1.395 1.966 



NACA AIRFOIL 12 

NUMBER OF PANELS USED = 20 




-8.0 -4.0 0.0 4.0 8.0 12.0 16.0 

ANGLE OF ATTACK 

Figure 18. Cl & Cm c/4 vs. Alpha 



50 



APPENDIX C 

PROGRAM NEW_VOR USER’S MANUAL 
USER’S GUIDE CONTENTS 

I. INTRODUCTION 52 

II. ASSUMPTIONS AND LIMITATIONS 52 

III. INPUT DESCRIPTION 52 

IV. INPUT RESTRICTIONS 53 

V. SAMPLE PROBLEMS 53 

VI. STARTING THE PROGRAM 53 

VII. SAMPLE GRAPHICAL OUTPUTS 54 



51 



I. INTRODUCTION 



The purpose of the NEW_VOR program is to provide an application of 
the vortex lattice method for the determination of the lift distribution of a flat 
rectangular plate. This method is based on a distribution of discrete horse- 
shoe vortices over a wing surface that has been divided into a finite number 
of panels. A system of linear equations is developed for the vortex strengths 
on the pemels and solved by matrix methods. 

II. ASSUMPTIONS AND LIMITATIONS 

This program is limited to flat rectangular wings which it divides into 
panels, using a uniform grid. Additionally, the uniform grid spacing method 
incorporates an enhancement whereby the panels do not extend to the wing 
tips, but only to a distance of d/4 from the tips. The value of d is the spanwise 
width of a wing panel. 

The solution is determined for conditions of incompressible and inviscid 
irrotational flow. Since we are considering an inviscid fluid, the coefficient of 
drag provided in the results is an accumulation of numerical errors. This pro- 
gram is intended to be used for the emalysis of flat rectangular wings with low 
aspect ratio. 

in. INPUT DESCRIPTION 

There are very few input values required for this simple program. Their 
description and program variable names are listed below. 

AR — ^Aspect ratio of the wing. (Span) 2/Area or Span/Chord. 

NXrNY — ^Number of vortices in the X and Y directions. 

ALPHA — ^Angle of attack. (Angle between the chord and the freestream velocity.) 
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IV. INPUT RESTRICTIONS 

The program, as written, is limited to 350 total surface vortices. This 
may be modified by changing the size of the arrays, however for the wings 
that this program was intended to analyze, this should be sufficient. The pro- 
gram will accept values for ALPHA up to 45 degrees, but, as noted previously 
with program NEW_PANEL, the user is cautioned that values above about 
10° may result in output data which in incorrect. 

V. SAMPLE PROBLEMS 

Two sample problems will be used to illustrate the use of the NEW_VOR 
program. The first nm will use a flat rectangular wing with an aspect ratio of 
two. The lattice will be created by placing three vortices on the wing in the X 
direction and five vortices on the wing in the Y direction. The vortices will be 
distributed using the Uniform Grid spacing method and the wing will be set 
at an angle of attack (alpha) of six degrees. The second run will use the same 
wing, but wdth five vortices on the wing in the X direction and 10 vortices on 
the wing in the Y direction. 

VI. STARTING THE PROGRAM 

Begin with the screen showing the DCL prompt, which looks like this: 

$ 

Next, ensure that the program is in your directory by typing: 

DIR [Return] 

and viewing the files for NEW_VOR.EXE. 

To nm the program, type: 

RUN NEW_VOR [Return] 
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VII. SAMPLE GRAPHICAL OUTPUTS 
A. SAMPLE PROBLEM ONE 

The program will start and the screen should look similar to what is 
shown in Figure 19. 



PROGRAM VORLAT : VERSION 5:10 OCTOBER 89 

VORTEX-LATTICE METHOD USED TO DETERMINE SPANWISE 
LIFT DISTRIBUTION FOR A FLAT RECTANGULAR WING 

ENTER THE ASPECT RATIO? 



Figure 19. Initial Screen for Program NEW_VOR 



Respond to the request for the aspect ratio by entering: 

2 [Return] 

Respond to the request for the number of vortices by entering: 

3,5 [Return] 

Finally, enter the angle of attack in degrees: 

6 [Return] 

The screen is then cleared and you will be presented with what is shown in 
Figure 20. If your display agrees with this, respond to the question by 
entering: 

1 [Return] 
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THE CURRENT VALUES ARE; 

1 ) ASPECT RATIO = 2.000000 

2) NUMBER OF VORTICES (NX, NY) = 3, 5 

3) ANGLE OF ATTACK (DEGREES) = 6.000000 

THE CALCULATED PARAMETERS ARE: 

DELTA X = 0.3333333 
DELTA Y = 0.1904762 

NUMBER OF EQUATIONS TO SOLVE = 15 
ARE THESE VALUES CORRECT? (YES=1, NO=2) 



Figure 20. Data Review/Correction Screen 



If you should desire to change any values, enter 2, and you will be 
asked which value you want to correct and the new desired value. Following 
entry of the correct values and a positive response, the program begins the 
solution process. It returns with the coefficients of lift and drag at the indi- 
cated spanwise positions, as well as the chordwise center of pressure for those 
positions. Overall values for the coefficients of lift, drag, induced drag and 
moment about the leading edge are calculated and then printed out near the 
bottom of the screen. Don't worry if you miss some of the values as they scroll 
up on the screen. All the values are printed to both the screen and to the data 
file (VORLAT4.DAT). 

The program now asks if you wemt to print the results. Entering 
£ui affirmative response of 'Y' will print the output file VORLAT4.DAT. 

The program will now ask if you want to graph the results. Enter: 

Y [Return] 
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Your screen should now look similar to Figure 21. 



WHICH OF THE FOLLOWING RELATIONSHIPS 
DO YOU WANT TO GRAPH? 

1) CLVS.Y 

2) CD VS.Y 

3) CL VS. CD 

4) NONE 

INPUT OPTION NO. (1,2,3 OR 4) 



Figure 21. Graphical Selection Screen 

Enter a desired plot selection and compare your one plot to the 
sample output plot at the end of this section. There should not be any differ- 
ence. You will also be asked if you would like a print of the respective plot. 
Upon entering: 

N [Return] 

your screen should once again be similar to Figure 21. Enter a response of 4 
to exit the graphing loop. 

B. SAMPLE PROBLEM TWO 

The program now asks if you want to make another run. Enter: 

1 [Return] 

You should now be back at the data review/correction screen and it should 
look like Figure 20. 
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Now run the same wing, but change the number of vortices to 5 
and 10. Enter: 

2 [Return] 

You want to change the number of vortices, so enter 

2 [Return] 

Respond to the request for the number of vortices by entering: 

5,10 [Return] 

The screen is automatically updated and you will see that the 
number of vortices has changed. As in the previous example, responding with 
a 'T causes the program to proceed to the output stage. The solution will be 
printed to the screen and appended to the data file which contains the data 
from the prior run. Again you will be afforded the opportunity to print and 
graph the results as in Sample Problem One. Respond accordingly... the out- 
put file/graphical plots for all plotting selections are enclosed at the end of 
this section. 

The program now asks if you want to make another run. The ses- 
sion is finished, so enter 

2 [Return] 

This completes the sample problems for the NEW_VOR program. 
Figures 22 through 24 give representative graphical outputs created by these 
sample problems. To create these plots, five vortices across the wing chord 
(NX) and 10 vortices across the span (NY) were used. 
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Y/B/2 

Figure 22. Cl vs. Y 
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2-D PLOT 
• = CD VALUES 



FLAT RECTANGULAR WING 
ASPECT RATIO (AR) = 2.00 

ANGLE OF ATTACK = 6. 00 

1 1 1 



^0. 000 0. 200 0. 400 0. 600 
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Figure 23. Cd vs. Y 
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DEGREES 

— I 

0. 800 1. 



p. 000 0. 065 0. 129 0. 194 0. 258 0. 323 




Figure 24. Cl vs. Cd 
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I. INTRODUCTION 



The SUB program has been adapted from a National Aeronautics and 
Space Administration (NASA) FORTRAN program and has been used consid- 
erably at the Langley Research Center. Additionally, this particular program 
has also been used in industry and the results have shown good correlation 
with experimental values. SUB has subsequently been revised to enhance its 
ease of use and its abihty to present accurate graphical results. 

The purpose of the SUB program is to estimate the subsonic aerody- 
namic characteristics of complex planforms. The program represents a lifting 
planform with a vortex lattice. A relatively complex planform may be analy- 
zed by creating the planform with up to 24 line segments on a semispan. 
Additionally, these line segments may have an outboard variable-sweep 
pamel or they may have several dihedral angles across the span. Furthermore, 
two planforms may be used together to represent a combination of wings and 
tails or wing, bodies, and tails. 

II. ASSUMPTIONS AND LIMITATIONS 

The use of this program is confined to the subsonic flow regime. Addi- 
tionally, the planform is in steady, uniform, inviscid, incompressible, attached 
flow conditions. 

Certain restrictions must also be kept in mind when using this program. 
Three specific restrictions apply to all planforms analyzed: 1) Only a total of 
two planforms may be specified; 2) The maximum number of horseshoe 
vortices on the left side must be limited to 120. When two planforms are 
specified, the sum total of the vortices is limited to 120. Within this limit, the 
number of horseshoe vortices in any chordwise row may vary from 1 to 20 and 
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the number of chordwise rows may vary from 1 to 50, and 3) The left side of 
the planform must be described with less than 24 line segments. 

Additionally, there are also three limitations which must be applied to 
variable-sweep planforms: 1) There should always be a fixed-sweep panel be- 
tween the root chord and the outboard variable-sweep panel; 2) The pivot 
cannot be canted from the vertical, and 3) Dihedral considerations cannot be 
programmed for the variable-sweep panel or at the intersection of this panel 
with the fixed portion of the wing. 

Finally, there exists three limitations when considering planforms which 
have nonzero dihedral angles or to two planforms which do not lie in the 
same plane: 1) The variation in local chord must be continuous from the tip 
chord to the root chord of each planform specified; 2) The number of horse- 
shoe vortices in each chordwise row must be at least two, and 3) The number 
of horseshoe vortices must be constant over the semispan of each planform. 

III. INPUT DESCRIPTION 

There are relatively few input values required for this program. Their 
description and program variable names are listed below. The user’s first task 
before ninning this program will be to create an input data file corresponding 
to the respective planform to be analyzed and the desired program specifica- 
tions. Each line of the input file is detailed explicitly. 
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A. GROUP ONE DATA 
Line 1 



PLAN Number of planforms for the configuration, 

TOT AL Number of sets of group two data (normally one). 

CREF Reference chord of the configuration (greater than zero). 

SREF Reference area of the configuration (greater than zero). 

Line 2 

AAN(IT) Number of line segments used to describe the left half of 
the planform. 

XS(IT) X location of the pivot; use 0 on a fixed wing, 

YS(IT) Y location of the pivot; use 0 on a fixed wing. 

RTCDHT(IT) Vertical distance of the particular planform being 

read in with respect to the wing root chord height; use 0. 

** The next series of input data lines are used to describe each line segment 
which was used to specify the planform shape. In other words, if one has used 
five line segments to describe his or her planform, the next five lines will des- 
cribe each line segment respectively. The first break-point is located at the 
intersection of the left wing leading edge with the root chord. They are num- 
bered in increasing order for each intersection of lines in a counterclockwise 
direction. The input variables for each of these lines is as follows; 

Lines 3-Whatever 

XREG(I,IT) X location of the ith breakpoint. 

YREG(I,IT) Y location of the ith breakpoint. 

DIH(I,IT) Dihedral angle (degrees) in y-z plane of line from 
breakpoint; positive upward. 
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AMCD The move code. (This input indicates whether or not the 
line segment in question is on a movable panel. Use 1 
for a line which is fixed or 2 for a line which is movable. 



B. GROUP TWO DATA 

Depending upon planform specifications, group two data could con- 
sist of three sections. The first section is always included. The second section 
is to be used if the number of chordwise horseshoe vortices varies across the 
semispan. One reason to vary the number of chordwise vortices across the 
semispan would be the need to analyze specific portions of the wing which 
may experience great pressure gradients i.e. at the intersection point of the 
fixed and the movable planforms for a variable sweep wing. The third section 
is used when the wing has twist and/or camber distribution and may consist 
of up to 15 lines, depending upon the number of horseshoe vortices. 



Line 1 (Section One) 

CONFIG An arbitrary configuration number (up to four digits) — 
user’s choice. 

sew The number of chordwise horseshoe vortices to be used to 

represent the wing; a maximum value of 20 may be used. 
If the user desires that the number of chordwise vortices 
vary across the semi span, enter 0. Entering zero will 
require the use of section two of Group Two data. The 
sew = 0 option can only be used on wings without dihe- 
dral and for coplanar wing-tail configurations. 

VIC The number of spanwise rows at which chordwise horse- 

shoe vortices will be TBLSCWd) cannot exceed 120. 

MACH Mach number. A value other than 0 will cause the 
Prandtl-Glauert compressibility factor to be applied. Re- 
gardless, the Mach number should be less than the critical 
Mach number. 
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CLDES Desired Lift Coefficient. Used to obtain the span load dis- 
tribution at a particular lift coefficient. If this aspect is 
not required enter 1. Enter 11 for drag polar data. 

PTEST If the damping-in-roll parameter is desired, enter 1. 

QTEST If CLq or Cmq stability derivatives are desired, enter 1. 

However, PTEST and QTEST cannot both be done in the 
same program nm. 



TWIST(l) Twist code for the first planform. Enter 0 for no twist. 

Enter 1 if the planform has twist and provide data in sec- 
tion three. 



SA(1) Variable sweep angle for the first planform. Specify the 

leading-edge sweep angle (degrees) for the first movable 
line adjacent to the fixed portion of the planform. For a 
fixed planform, this quantity may be omitted. 

TWIST(2) Twist code for the second planform.** 

SA(2) Variable sweep angle for the second planform.** 



**Obviously, these inputs may be omitted if there is only one planform. 



Line 2 (Section Two) 

Again, section two is to be used if SCW was set to 0 thus allowing for the 
number of chordwise horseshoe vortices to vary across the semispan. 

STA Total number of spanwise rows of horseshoe vortices per 

semispan. This input sets the number of values to be read 
into TBLSCWd) — next input. 

Lines 3- Whatever 

TBLSCWd) Number of horseshoe vortices in each row starting at 
the row near the tip of the first planform and proceeding 
to the row near the root. If a second planform has been 
specified, the table of chordwise rows concludes with the 
number of vortices specified for the second planform (see 
Example B for format). 
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Section Three 

Again, section three is to be used if the planform has twist. 

ALP(NV) Local angle of attack in radians. Refer to Example Three 
(3) 

FORMAT Refer to the sample input data files on how to properly 
format the input data files. Failure to follow these exam- 
ples implicitly will result in a "data read error". 

IV. SAMPLE PROBLEMS 

Three sample problems have been included in this user's guide section. 
The first two problems analyze a fixed planform without a variable-sweep 
panel. The first problem simply uses four (4) spanwise vortices while the sec- 
ond problem uses 40 spanwise vortices. The second problem demonstrates the 
benefit of using extra horseshoe vortices to enhance data representation. The 
last problem is an example of a rather complex planform. This particular 
wing has variable chordwise vortices across seven (7) spanwise rows and has 
twist incorporated into the wing. Additionally, this wing is described using 14 
line segments. 

V. STARTING THE PROGRAM 

Begin with the screen showing the DCL prompt, which looks like this: 

$ 

Next, enter the following command: 

SET DEF [.SUB] 

Now, enter the command to run the program: 

RUN SUB 
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The program will start £ind the screen should look similar to what is shown is 
Figure 25. 



PRCX3RAM MLVL - SUBSONIC VORTEX LATTICE ANALYSIS 
ENTER INPUT DATA FILE NAME 

USE LAST.END AS DATA FILE NAME TO STOP THE PROGRAM 



Figure 25. Initial Screen for Program SUB 



VI. SAMPLE GRAPHIC OUTPUTS 
A. EXAMPLE PROBLEM 1 

Enter the name of the input data file. 

A9WS60.DAT [Return] 

Once the program has finished its data tabulations, your screen 
should be similar to Figure 26. 



PROGRAM RESULTS HAVE BEEN WRITTEN TO THE FILE 
OUTFILE.DAT. 

WOULD YOU LIKE A PRINTED COPY OF THIS OUTPUT FILE? 
YES OR NO (Y/N) 



Figure 26. Printing Determination Screen 
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Respond negatively to this request by typing: 

N [Return] 

Respond affirmatively to the request to copy the output data file 
(OUTFELE.DAT) to another file by typing: 

Y [Return] 

A screen similar to Figure 27 will then appear which lists the file 
choices possible for copying. 



WHAT NAME WOULD YOU LIKE FOR THE OUTPUT FILE? 


1) 


VIGILANTE.DAT 


2) 


CORSAIR.DAT 


3) 


HAWKEYE.DAT 


4) 


SKYHAWK.DAT 



Figure 27. Output File Designation Screen 



Select from the designated list of file names your choice. 

1,2,3, OR 4 [Return] 

Respond affirmatively to the request to graph the results by typing: 

Y [Return] 

The screen should now look like Figure 28. 
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WHICH OF THE FOLLOWING RELATIONSHIPS 
DO YOU WANT PLOTTED? 


1) 


INDUCED DRAG COEFF VS. 2Y/B 


2) 


LE EDGE THRUST COEFF VS. 2Y/B 


3) 


SUCTION COEFF VS. 2Y/B 


4) 


SPAN LOAD COEFF VS. 2Y/B 


5) 


CL RATIO VS 2Y/B 


6) 


NONE 

INPUT OPTION NO. (1 .2.3, 4, 5, OR 6) 



Figure 28. Plot Determination Screen 



Select from the designated list of graphical relationships your 

choice. 

1,2, 3, 4, OR 5 [Return] 

The requested plot will momentarily appear on your screen. If you 
have remoted your terminal to terminal "KELLY" for printing purposes your 
plot will come up on the "ILELLY" monitor. Compare your plot with the 
example plots corresponding to EXAMPLE 1; it should be the same. 

Respond negatively to the request to print the plot by typing: 

N [Return] 
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The user will again be given the opportunity to graph another rela- 
tionship (Figure 28 wiU be presented). Respond with the 6th choice to exit the 
graphing loop. Enter: 

6 [Return] 

The program now asks if you want to make another nm. Enter 

1 [Return] 

B. EXAMPLE PROBLEM 2 

The screen should again look like Figure 25. 

Enter the name of the input data file. 

E9WS60.DAT [Return] 

Respond negatively to the request for a printed copy of the output 
file by typing: 

N [Return] 

Respond negatively to the request to copy the output data file 
(OUTFILE.DAT) to another file by typing: 

N [Return] 

Respond affirmatively to the request to graph the results by typing: 

Y [Return] 

Again, Figure 28 will appear on your screen with a listing of the 
available plotting routines. 

Select from the list yoiu- plotting choice. 

1,2, 3, 4, OR 5 [Return] 
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The requested plot will momentarily appear on your screen. Again, 
if you have remoted your terminal to terminal "KELLY” for printing purposes 
your plot will come up on the "KELLY' monitor. Compare your plot with the 
example plots corresponding to EXAMPLE 2; it should be the same. 

Respond negatively to the request to print the plot by typing: 

N [Return] 

The user will again be given the opportunity to graph another rela- 
tionship. Respond with the 6th choice to exit the graphing loop. Enter: 

6 [Return] 

Respond affirmatively to the request to perform another run of pro- 
gram by typing: 

1 [Return] ** 

**Entering a “2” would exit the user from the program. 

C. EXAMPLE PROBLEM 3 

Enter the name of the input data file. 

B9WS60.DAT [Return] 

Respond negatively to the program request to print OUTFILE.DAT. 

Enter: 

N [Return] 

Respond negatively to the request to copy the output data file 
(OUTFILE.DAT) to another file by typing: 

N [Return] 
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Respond affirmatively to the request to graph the results by typing: 

Y [Return] 

Select from the designated list of graphical relationships your 

choice. 

1,2,3, 4, OR 5 [Return] 

The requested plot will then appear on your screen and you will be 
asked if you want to print the plot. Compare your plot with the example plots 
corresponding to EXAMPLE 3; it should be the same. 

Respond negatively to the request to print the plot by typing: 

N [Return] 

The user will again be given the opportunity to graph another rela- 
tionship. Respond with the 6th choice to exit the graphing loop. Enter: 

6 [Return] 

Respond negatively to the request to perform another run of pro- 
gram by typing: 

2 [Return] 

This completes the sample problems for the SUB program. Graphi- 
cal output examples created by these sample nms are shown in Figures 29 
through 34. The first five plots were generated from the analysis of a wing 
with an aspect ratio of nine and a leading edge sweep angle of 60°. The last 
plot (Figure 34) was produced from the analysis of a rather complex planform 
[Ref. 6], which had seven rows of spanwise vortices with nine vortices across 
the chord at horseshoe vortex Number 3. 
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INDUCED DRAG COEFFICIENT 

-0.021 -0.010 0.001 0.012 0.022 0.033 



CONTRIBUTION TO TOTAL COEFF. 
INDUCED DRAG COEFFICENT 




-| 1 1 1 1 1 

.000 0.200 0.-100 0.600 0.800 1.000 

2Y/B 



Figure 29. Induced Drag Coeff. vs. 2Y/B 
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LE THRUST COEFFICIENT 

-0.015 -0.006 0.004 0.013 0.023 0.032 0.042 




. 000 0. 200 0. 400 0. 600 0. 800 1. 000 

2Y/B 



Figure 30. LE Thrust Coeff. vs. 2Y/B 
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SUCTION COEFFICIENT 

- 0.020 - 0.001 0.018 0.037 0.056 0.075 0.094 0.113 




. 000 0 . 200 0 . 400 0 . 600 0 . 800 1 . 000 

2Y/B 



Figure 31. Suction Coeff. vs. 2Y/B 



76 



SPAN LOAD COEFFICIENT 

0.531 0.641 0.751 0.861 0.971 1.081 




2Y/B 



Figure 32. Span Load Coeff. vs. 2Y/B 
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COEFFICIENT OF LIFT RATIO (SECTION/WING) 




2Y/B 



Figure 33. Coeff. of Lift Ratio vs. 2Y/B 
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COEFFICIENT OF PRESSURE CHANGE 

•0. 420 0. 972 1. 525 2. 077 2. 629 3. 182 



HORSESHOE VORTEX 



NUMBER 3 




Figure 34. Delta Cp vs. X c/4 
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I. INTRODUCTION 



The SUPER program has been adapted from a National Aeronautics and 
Space Administration (NASA) FORTRAN program and has been used consid- 
erably at the Langley Research Center. Additionally, this particular program 
has also been used in industry. The results have shown good correlation with 
experimental results. SUPER has subsequently been revised to enhance it's 
ease of use and its ability to present accurate graphical results. 

The purpose of the SUPER program is to estimate the supersonic aero- 
dynamic characteristics of complex planforms. Linearized supersonic lifting 
surface theory is employed to calculate the aerodynamic characteristics of a 
warped wing of arbitrary planform. The program calculates lifting pressure 
distribution for the chordwise warped wing at fixed attitude and the pressure 
distribution (per degree angle of attack) for a corresponding flat wing. These 
two pressure distributions are combined by superposition principles and inte- 
grated over the wing surface to obtain the variation of aerodynamic charac- 
teristics with changes in angle of attack. 

II. ASSUMPTIONS AND LIMITATIONS 

The use of this program is confined to the supersonic flow regime. In ad- 
dition, the Hnearized supersonic lifting surface theory, used in this program, 
applies to wings having negligible thickness. 

There exist two specific limitations which must be considered when en- 
tering the respective input data values. The number of semispan grid ele- 
ments is limited to 100 or 47.5*B*SPAN/XMAX. The relative increase in 
semispan grid elements will increase the computational time of the program. 
Additionally, the number of percent chord values is limited to 26. Lastly, 
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there are a few other input restrictions which need to be referenced when 
creating your input data file. The next section delineates each respective 
input and declares auiy restrictions. 

III. INPUT DESCRIPTION 

There are relatively few input values required for this program. Their 
description and program variable names are listed below. The user's first task 
before running this program will be to create an input data file corresponding 
to the respective planform to be analyzed and the desired program 
specifications. 

LINE1: $INPT1 

Type line as indicated. This lines cues the program for input of data. 

LINE 2: XM 

Mach Number of freestream. 

LINE 3: NOM 

Number of additional Mach Numbers other than XM (NOM<5). 

LINE 4: NOPCT 

Number of percent chord values for TZORD input (NOPCT<26). 

LINE 5: TPCT 

Table of percent chord values, corresponding to NOPCT, in increasing 
order from 0 to 100. 

LINE 6: JBYMAX > 

Number of spanwise stations at which TZORD is to be specified (JBYMAX<51). 
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LINE 7: 


TYB2 

Table of semispan fractions, corresponding to JBYMAX, in increasing order from 
0 to 1 .0. 


LINE 8: 


TZORD 

Zq coordinates of right-hand wing panel corresponding to TYB2 and TPCT. All values 
of at a given semispan station entered in order according to TPCT, 26 values 
required to fill a table column. You must enter 26 values per column although only 
NOPCT values are used. After the first column is filled, repeat with other TYB2 
stations, proceeding to right-hand wing tip. 


LINE 9: 


REFAR 

Wing reference area. 


LINE 10: 


SPAN 

Total wing span. 


LINE 11: 


XLEO 

X coordinate of wing leading edge of y=0. 


LINE 12: 


XTEO 

X coordinate of wing trailing edge at y=0. 


LINE 13: 


XMAX 

Largest value of x in wing definition. 


LINE 14: 


xo 

Distance from some arbitrary location to wing apex. XO=0. if you are considering the 
wing only. This term is used in locating streamwise lift distribution with respect to XO 
rather than wing apex. 


LINE 15: 


TYPEX 

= 0. Input TTXLE and TXTE tables. 

= 1. Input NLEX, NTEX and tables of TBLEX, TBLEY, TBTEX, TBTEY. 
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LINE 16: TXLE 

Table of wing leading edge x coordinates at successive values of 
y=((SPAN/2)/NON)*N where N = 1 , 2. 3, ...NON. (Omit if TYPEX = 1 .) 

LINE17:TXTE 

Table of wing trailing edge x coordinates specified at same values of y as TXLE. 
(Omit if TYPEX = 1.) 

LINE 18: NLEX 

Number of leading edge (x;y) points to be input (NLEX<15). (Omit if TYPEX = 0.) 

LINE 19: NTEX 

Number of trailing edge (x,y) points to be input (NTEX<15). (Omit if TYPEX = 0.) 

LINE 20: TBLEX 

Table of NLEX leading edge x values(spanwise, root to tip). (Omit if TYPEX =0.) 

LINE 21: TBLEY 

Table of NLEX leading edge y values(spanwise, root to tip). (Omit if TYPEX =0.) 

LINE 22: TBTEX 

Table of NTEX trailing edge x values(spanwise, root to tip). (Omit if TYPEX =0.) 
LINE 23: TBTEY 

Table of NTEX trailing edge y values(spanwise, root to tip). (Omit if TYPEX =0.) 

LINE 24: CBAR 

Reference length used for pitching moment coefficient. 

LINE 25: XMREF 

X distance from X=0. locating pitching moment center. 
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LINE 26: NON 

Number of semispan grid elements selected to represent the wing. (NON^O or 
NON^7.5‘B*SPAN/XMAX (whichever value is less). 

LINE 27: $END 

Line statenrvent ends input of data. 

IV. SAMPLE PROBLEMS 

Two sample problems have been included in this user's guide section. 
Both consider the same planform shape, but the input method of the plan- 
form shape is different. Only one set of plots exists in the sample problems 
output file section in that the two sets of plots are identical. 

V. STARTING THE PROGRAM 

Begin with the screen showing the DCL prompt, which looks like this: 

$ 

Next, enter the following command: 

SET DEF [.SUPER] [Return] 

Now, enter the command to run the program: 

RUN SUPER [Return] 

The program will start and the screen should look similar to what is 
shown is Figure 35. 
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PROGRAM A4410 - SUPERSONIC VORTEX LATTICE ANALYSIS 



ENTER THE INPUT FILE NAME 
USE LAST.END AS THE DATA FILE NAME 
TO STOP THE PROGRAM 



Figure 35. Initial Screen for Program SUPER 



VI. SAMPLE GRAPHICAL OUTPUTS 
A. EXAMPLE PROBLEM 1 

Enter the name of the input data file. 

SSVL1.DAT [Return] 



Once the program has finished its data tabulations, your screen should be 
similar to Figure 36. 



PROGRAM RESULTS HAVE BEEN WRITTEN TO THE FILE 
OUTFILE.DAT. 

WOULD YOU LIKE A PRINTED COPY OF THIS OUTPUT FILE? 
YES OR NO (Y/N) 



Figure 36. Printing Determination Screen 



Respond negatively to print request by typing; 

N [Return] 
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Respond affirmatively to the request to copy the output data file 
(OUTFILE.DAT) to another file by typing: 

Y [Return] 

A screen similar to Figure 37 will then appear which lists the file choices pos- 
sible for copying. 



WHAT NAME WOULD YOU LIKE FOR THE OUTPUT FILE? 

1) TOMCAT.DAT 

2) PHANTOM.DAT 

3) INTRUDER.DAT 

4) CRUSADOR.DAT 

ENTER 1 , 2, 3 OR 4 



Figure 37. Output File Designation Screen 



Select from the designated list of file names your choice. 

1,2,3, OR 4 [Return] 

Respond affirmatively to the request to graph the results by typing: 

Y [Return] 

The screen should now look like Figure 38. 
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WHICH OF THE FOLLOWING RELATIONSHIPS 
DO YOU WANT PLOTTED? 



1) SPANWISE PRESSURE DISTRIBUTION 

2) CHORDWISE PRESSURE DISTRIBUTION 

3) DRAG POLAR (CL VS. CD) 

4) STREAMWISE LIFT DISTRIBUTION 

5) SPANWISE LIFT DISTRIBUTION 

6) NONE 

INPUT OPTION NO. (1, 2, 3, 4, 5 OR 6) 



Figure 38. Plot Determination Screen 



Select from the designated list of graphical relationships your choice. 

1,2, 3, 4, OR 5 [Return] 

The requested plot will then appear on your screen and you will be 
asked if you want to print the plot. If you have remoted your terminal to 
terminal "KELLY" for printing purposes your plot will come up on the 
"KELLY' monitor. Compare your plot with the example plots corresponding 
to SAMPLE # 1; it should be the same. 

Respond negatively to the request to print the plot by typing; 

N [Return] 
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The user will again be given the opportunity to graph another rela- 
tionship (Figure 22 will be presented). Respond with the 6th choice to exit the 
graphing loop. Enter: 

6 [Return] 

The program now asks if you want to make another run. Enter 

1 [Return] 

B. EXAMPLE PROBLEM 2 

The screen should again look like Figure 35. 

Enter the name of the input data file. 

EXPROB2.DAT [Return] 

Respond negatively to the request for a printed copy of the output 
file by typing: 

N [Return] 

Respond negatively to the request to copy the output data file 
(OUTFILE.DAT) to another file by typing: 

N [Return] 

Respond affirmatively to the request to graph the results by typing: 

Y [Return] 

Again, Figure 38 will appear on your screen with a listing of the available 
plotting routines. 
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Select from the list your plotting choice. 

1,2,3, 4, OR 5 [Return] 

The requested plot will appear on your screen. Again, if you have 
remoted your terminal to terminal "KELLY' for printing purposes your plot 
will come up on the "KELLY" monitor. Compare your plot with the example 
plots corresponding to SAMPLE#2; it should be the same. 

Respond negatively to the request to print the plot by typing: 

N [Return] 

The user will again be given the opportunity to graph another relationship. 
Respond with the 6th choice to exit the graphing loop. Enter: 

6 [Return] 

Respond negatively to the request to perform another run of program by 
typing: 

2 [Return] 

This completes the sample problems for the SUPER program. Graphical 
output examples created by these sample runs are shown in Figures 39 
through 42. These plots were created from the analysis of a B2 Bomber 
planform at a Mach of 1.2. The span used was 200 feet with a planform 
reference area of 8260.4 ft^. Thirty semispan grid elements were used to 
represent the wing. 
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COEFFICIENT OF PRESSURE 

0. 000 0. 008 0. 015 0. 023 0. 030 0. 038 0. 046 0. 053 



CHORDAL FRACTION (X/L) 



0. 690 




Figure 39. Spanwise Cp Distribution 
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COEFFICIENT OF PRESSURE 



SPAN FRACTION (Y/B/2) 



0. 398 




(X-XLE) /C 



Figure 40. Chordwise Cp Distribution 
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p. 000 0. 032 0. 063 0. 095 0. 127 0. 159 



DRAG POLAR CURVES 
• = DRAG POLAR - FLAT W/HG 
o = DRAG POLAR - CAMBERED WING 



0 . 000 



0. 100 




0. 400 



Figure 41. Drag Polar 



0. 500 
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LIFT FRACTION 

0. 000 0. 009 0. 017 0. 026 0. 035 0. 043 




0.000 0.200 0.400 0.600 

Y/B/2 



0. 800 



1. 000 



Figure 42. Spanwise Lift Distribution 
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APPENDIX F. PROGRAM NEW DOUBLE COMPUTER CODE 



APPENDIX F. PROGRAM NEW_DOUBLE COMPUTER CODE 
PROGRAM NEW_DOUBLE 

*** MODIFIED FOR USE ON THE MICROVAX/2000 BY J. A. CAMPBELL (JUL 88) 
UPDATES MADE BY C. M. MACALLISTER JAN-JUL 89 (CMM) 



INCOMPRESSIBLE AERODYNAMICS OF SYMMETRIC AIRFOIL 
AT ZERO ANGLE OF ATTACK BY LINE DOUBLET DISTRIBUTION 



ORIGINAL IBM MAINFRAME PROGRAM WAS ADAPTED FROM JACK MORAN'S BOOK 
'AN INTRODUCTION TO THEORETICAL AND COMPUTATIONAL AERODYNAMICS', 
WILEY AND SONS, NEW YORK 1984. THE LISTING IS FOUND ON PAGE 75. 



PROGRAM FLEXIBILITY AND USER INTERFACE WAS REVISED FOR 
PROFESSOR J. V. HEALEY BY JOHN CAMPBELL. 

ADDITIONAL PROGRAM UPDATES TO INCLUDE DUBLET USE FOR ANY 
ARBITRARY 2-D SHAPE, PRINTING ROUTINES, PROCESSING CORRECTIONS, 
AND GRAPHICAL ANALYSIS WERE MADE BY CRAIG MACALLISTER IN 
JAN-JUL 1989. (CMM) 






CHARACTER* 1 I ANS , PRINT , GRAPH , PLOT 1 , PLOT2 , 
+PLOT3 , CHECK , CORRECT 
INTEGER NANS , DATPO , PRINTOPT , GRAPHOPT 
REAL*4 T(100),M(100),XS,XF 
REAL XX, CP 
INTEGER N,R,NPRINT 
COMMON /GRAPH/XX,CP,NPRINT 
COMMON /MAIN/ T,M,N,XS,XF 
COMMON /GRAPHER/GRAPHOPT,XMAXY 
COMMON /FCN/AX,TAU,NTYPE 
COMMON /DATA/COORXC 101) ,COORY( 101) 

COMMON /PROB/DATPO 
DIMENSION NUM(IOO) 

REAL MPLOT 



OPEN FILE FOR DOUBLET STRENGTH DISTRIBUTION OUTPUT 
OPEN (UNIT=11, 

2 FILE= 'DUBLET.DAT', 

2 ORGANIZATION= 'SEQUENTIAL', 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE', 

2 FORM= 'FORMATTED', 

2 STATUS= 'UNKNOWN') 

OPEN FILE FOR BODY SHAPE OUTPUT 
OPEN (UNIT=12, 

2 FILE= 'SHAPE.DAT', 

2 ORGANIZATION= 'SEQUENTIAL' , 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE', 
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c 

c 



c 

c 



c 

c 



c 

c 

c 



2 FORM= 'FORMATTED', 

2 STATUS= 'UNKNOWN') 

OPEN FILE FOR BODY SURFACE PRESSURE DISTRIBUTION OUTPUT 
OPEN (UNIT=13, 

2 FILE= ' PRESSURE. DAT' , 

2 ORGANIZATION= 'SEQUENTIAL', 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE', 

2 FORM= 'FORMATTED', 

2 STATUS= 'UNKNOWN') 

OPEN ANOTHER FILE FOR BODY SURFACE PRESSURE DISTRIBUTION OUTPUT 
OPEN (UNIT=14, 

2 FILE= 'PRESS.DAT', 

2 ORGANIZATION= 'SEQUENTIAL', 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE', 

2 FORM= 'FORMATTED', 

2 STATUS= 'UNKNOWN') 

OPEN ANOTHER FILE FOR BODY SHAPE OUTPUT 
OPEN (UNIT=15, 

2 FILE= 'SHAPEBODY.DAT' , 

2 ORGANIZATION= 'SEQUENTIAL' , 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE', 

2 FORM= 'FORMATTED', 

2 STATUS= 'UNKNOWN') 



CALL LIBRARY ROUTINE TO CLEAR THE SCREEN, THE PRINT HEADER 
5 CONTINUE 
CALL CLRSCRN 



PRINT * 

PRINT *, ' PROGRAM DUBLET : VERSION 3 ; 4 OCTOBER 89 ' 
PRINT * 



* 

* 



PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT 
PRINT 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT 
PRINT 
PRINT 
PRINT 
READ (5 
IF (NTYPE 
PRINT 



* 

* 

* 

* 



*) 



DOUBLET DISTRIBUITON METHOD IS USED TO DETERMINE' 
INCOMPRESSIBLE AERODYNAMICS OF AN ELLIPSE, SYMMETRICAL' 
AIRFOIL OR ARBITRARY SYMMETRIC SHAPE AT ZERO ANGLE' 

OF ATTACK' 

f 

PROGRAM ASSUMES A NONDIMENS lONAL CHORD, THAT IS,' 

THE VALID RANGE OF X IS FROM 0 TO 1. ' 

ENTER TYPE OF BODY SHAPE DESIRED; ' 

1 ) E LL I PT I C ^ 

2) SYMMETRICAL AIRFOIL-LIKE OR' 

3) ARBITRARY SYMMETRIC SHAPE' 

ENTER 1, 2, OR 3. ' 

f 

NOTE THAT OPTION 3 WILL REQUIRE MANUALLY INPUTTING DATA' 
POINTS FOR THE UPPER SIDE OF THE RESPECTIVE BODY' 

NTYPE 

.LT. 1 .OR. NTYPE . GT. 3) THEN 

INVALID ENTRY. ENTER 1, 2, OR 3.' 



* ' 
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GO TO 15 
END IF 

IF (NTYPE .EQ. 3) THEN 
CALL CLRSCRN 

PRINT *, 'HOW MANY UPPER PROFILE DATA POINTS DO' 

PRINT 'YOU DESIRE? (ENTER A NUMBER BETWEEN 3 AND 100)' 
PRINT ' 

PRINT *, 'BE AWARE THAT THE LEADING EDGE OF YOUR DESIRED' 
PRINT *, 'SHAPE HAS BEEN PROGRAMMED TO BE AT THE ORIGIN' 
PRINT *, 'AND THAT YOUR TRAILING EDGE IS AT (1,0). SCALE' 
PRINT *, 'YOUR SHAPE/OBJECT ACCORDINGLY. ' 

PRINT ' 

17 READ (5,*) DATPO 

IF (DATPO .LT. 3 .OR. DATPO . GT. 100) THEN 

PRINT *, 'INVALID ENTRY. ENTER A NUMBER BETWEEN' 

PRINT *, 'THREE(3) AND 100 INCLUSIVE. ' 

GO TO 17 
END IF 

DO 26 R = 1, DATPO 
COORX(l) =0.0 
COORX(DATPO+2) =1.0 
WRITE (5,27) R 

27 FORMAT (IX, 'ENTER X(',I2,')') 

READ (5,*) COORX(R+l) 

COORY(l) = 0. 0 
COORY(DATPO+2) = 0. 0 
WRITE (5,28)R 

28 FORMAT (IX, 'ENTER Y(',I2,')') 

READ (5,*) COORY(R+l) 

26 CONTINUE 

PRINT *, ' ' 

PRINT *, ' WOULD YOU LIKE TO CHECK YOUR SURFACE DATA POINTS? ' 

PRINT*, ' (Y/N)' 

READ 1002, CHECK 

IF (CHECK. EQ. 'Y' . OR. CHECK. EQ. 'y' )THEN 
313 CALL CLRSCRN 

DO 65 I = l,DATPO+2 

WRITE(5,29) I,COORX(I) ,COORY(I) 

29 FORMAT(5X,I3,3X,F8.4,3X,F8. 4,/) 

65 CONTINUE 

PRINT *, ' WOULD YOU LIKE TO MAKE ANY CORRECTIONS?' 

PRINT *, ' (Y/N)' 

PRINT *, ' ' 

READ 1002, CORRECT 

IF (CORRECT. EQ. 'Y' . OR. CORRECT. EQ. 'y' )THEN 
PRINT *, ' ' 

PRINT *, ' WHICH DATA POINT WOULD YOU LIKE TO CORRECT?' 
NUMBERS = DATPO + 2 
WRITE (5, 30) NUMBERS 

30 FORMAT (5X, 'ENTER A NUMBER 1 THRU', 14,' INCLUSIVE') 

312 READ (5,*)NUMCOR 

IF (NUMCOR. LT. 1. OR. NUMCOR. GT. NUMBERS)THEN 

PRINT *, ' INVALID ENTRY ' 

WRITE (5,30)NUMBERS 
PRINT *, ' ' 

GO TO 312 
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ENDIF 

WRITE (5,27)NUMCOR 
READ ( 5 , * ) COORX( NUMCOR ) 

WRITE (5, 28) NUMCOR 
READ(5,*)C00RY(NUMC0R) 

GO TO 313 
ENDIF 
ENDIF 
GO TO 70 
END IF 

PRINT ENTER THICKNESS RATIO (TAU). ' 

READ (5,*) TAU 
IF (NTYPE .GT. 1) THEN 
PRINT * 

PRINT ENTER THE NONDIMENSIONAL X LOCATION OF MAXIMUM' , 

+ ' THICKNESS. ' 

20 READ (5,*) XMAXY 

IF (XMAXY .GT. 0.5) THEN 

PRINT THE PROGRAM CONSIDERS THE ONSET FLOW TO BE* 

PRINT APPROACHING FROM THE LEFT. THEREFORE, THE' 

PRINT X LOCATION OF MAXIMUM THICKNESS MUST BE < 0. 5. ' 

PRINT ~> PLEASE REENTER. ' 

GO TO 20 
END IF 

AX = (.5 * TAU)/(SQRT(XMAXY)*(1. - XMAXY)) 

END IF 

INPUT NUMBER OF INTERVALS N 

70 CALL CLRSCRN 
PRINT * 

PRINT *, ' ENTER NUMBER OF INTERVALS DESIRED. N =' 

71 READ (5,*) N 
PRINT * 

IF(N .LT. 2 .OR. N . GT. 100) THEN 
WRITE(6,21) N 

PRINT ' A MINIMUM OF TWO INTERVALS AND A MAXIMUM OF' 

PRINT ' 100 IS ALLOWED. ===> PLEASE REENTER. ' 

GO TO 71 
END IF 

21 FORMATC IX, 5X, 'NUMBER OF INTERVALS REQUESTED =',I3) 



ASK USER FOR AUTOMATIC OR MANUAL DETERMINATION OF ENDPOINTS. 
80 CONTINUE 
CALL CLRSCRN 
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PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
READ (5 



IF (NMETH 



WHICH METHOD DO YOU WISH TO USE TO DETERMINE THE' 
DOUBLET DISTRIBUTION ENDPOINTS? ' 

1) PROGRAM INTERVAL HALVING SUBROUTINE TO ITERATE. 

2) MANUAL ITERATION BY THE USER. ' 

3) RETURN TO START' 

4) EXIT PROGRAM' 

ENTER 1,2,3 OR 4' 



*) 



NMETH 
LT. 1 



OR. NMETH .GT. 4) THEN 



! 
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PRINT * , ' ' 

PRINT *, ' INVALID ENTRY. ENTER A NUMBER BETWEEN* 

PRINT *, ' ONE(l) AND F0UR(4) INCLUSIVE. ' 

GO TO 24 
END IF 

GO TO (120,100,5,999) NMETH 

MANUALLY DETERMINE ENDPOINTS OF SOURCE DISTRIBUTION XS , XF 

100 CONTINUE 

CALL CLRSCRN 
PRINT * 

PRINT *, ' ROUTINE FOR MANUAL DETERMINATION OF ENDPOINTS' 

PRINT * 

PRINT *, ' ' 

PRINT * 

PRINT *, ' ENTER THE DOUBLET DISTRIBUTION STARTING POINT, XS. ' 

PRINT *, ' (XS SHOULD BE APPROXIMATELY ONE HALF OF* 

PRINT *, ' THE NONDIMENS lONAL LEADING EDGE RADIUS.)' 

READ (5,*) XS 
PRINT * 

PRINT *, ' ENTER THE DOUBLET DISTRIBUTION ENDING POINT, XF. ' 

PRINT *, ' (XF SHOULD BE APPROXIMATELY ONE MINUS HALF* 

PRINT *, ' OF THE NONDIMENSIONAL TRAILING EDGE RADIUS. )' 

READ (5,*) XF 
PRINT * 

PRINT * 

CALL FINDM (T,M,N,XS,XF) 

CALL PRESS(0. 0,U0,CP0) 

CALL PRESS(1. 0,U1,CP1) 

GO TO 150 



120 CONTINUE 

CALL CLRSCRN 
PRINT * 

PRINT *, ' 

PRINT ' 
PRINT * 

PRINT *, ' 

PRINT * 



INTERVAL HALVING ROUTINE FOR DETERMINATION OF* 
DOUBLET DISTRIBUTION ENDPOINTS' 



t 



ENTER 
WHICH 
PRINT * 
PRINT * 
READ (5 
PRINT * 
PRINT * 



& 



THE PARAMETERS REQUIRED BY THE INTERVAL HALVING METHOD 
IS USED TO OBTAIN THE PROPER LOCATIONS FOR XS AND XF. 

' ENTER THE INTEGER EXPONENT FOR THE X TOLERANCE, NXTOL. 
' EXAMPLE: A VALUE OF 4, GIVES A TOLERANCE OF 0.0001. 

*) NXTOL 

' ENTER THE INTEGER EXPONENT FOR THE FUNCTION ' , 
'TOLERANCE, NFTOL. ' 



PRINT * 

READ (5 
PRINT * 

PRINT * 

PRINT * 

READ (5 
PRINT * 

PRINT *, ' ENTER THE OUTPUT PARAMETER, lOUT. ' 



' (SAME IDEA AS NXTOL; 5 YIELDS FTOL = 0.00001). ' 

*) NFTOL 

' ENTER THE MAXIMUM NUMBER OF ITERATIONS, MAXIT, TO ' 
' LOCATE XS AND XF, (FOR NFTOL = 6, SUGGEST 35-40) ' 
*) MAXIT 



I 

t 
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non 



c 



c 



PRINT *, ’ lOUT = 0 TO SUPPRESS ALL ITERATION RELATED OUTPUT' 

PRINT ’ 1 TO OUTPUT FINAL RESULTS ONLY' 

PRINT ' 2 TO OUTPUT DETAILS FOR EACH ITERATION' 

READ (5,*) lOUT 

CALL INTHV (NXTOL,NFTOL,NTYPE,MAXIT, I0UT,U0,U1) 

RUN THROUGH PROCESS AGAIN WITH FINAL VALUES OBTAINED BY ITERATION 
CALL FINDM (T,M,N,XS,XF) 

CALL PRESS(0. 0,U0,CP0) 

CALL PRESS(1. 0,U1,CP1) 



U AT X = 0 =' ,U0,' 
U AT X = 1 =' ,U1,' 



' THESE VALUES FOR U SHOULD BE NEAR ZERO. ' 
' DO YOU ACCEPT THESE RESULTS (Y/N)’ 



150 PRINT 
PRINT *, 

PRINT * 

PRINT 
PRINT *, 

READ 1000, IANS 
IF (IANS .EQ. 'N') THEN 

PRINT *, 'CORRECTION LINE NO. l' 
GO TO (120,100) NMETH 
ELSE 

GO TO 152 
END IF 



XS =' ,XS 
XF =' ,XF 



OUTPUT RESULTS 

152 PRINT 1010 

WRITE (11,1012) 

M(N+1) = 0. 0 

DO 200 I = 1,N+1 
MPLOT = REAL(M(I)*3. 1415926585) 
PRINT 1040, T( I), MPLOT 
200 WRITE (11,1040) T( I), MPLOT 
CLOSE (UNIT=11) 

PRINT 1020 
WRITE (12,1020) 

IF (NTYPE .LE. 2) THEN 
DO 210 I = 1,N 
XX = . 5*(T(I) + T(I+1)) 

YY = Y(XX) 

PRINT 1040, XX, YY 
WRITE (15,1040) XX, YY 

210 WRITE (12,1040) XX, YY 
XX = 1.0 

YY = 0. 0 

WRITE (15,1040) XX, YY 
END IF 

IF (NTYPE .EQ. 3) THEN 
DO 211 I = l,DATP0+2 
XX = COORX(I) 

YY = COORY(I) 

PRINT 1040, XX, YY 
WRITE (15,1040) XX, YY 

211 WRITE (12,1040) XX, YY 
END IF 

CLOSE (UNIT=12) 

CLOSE (UNIT=15) 
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PRINT 1030 

212 READ (5,*) NPRINT 

IF (NPRINT . LT. 2) THEN 

PRINT ' YOU MUST ENTER A MINIMUM OF 2. PLEASE REENTER. ' 
GO TO 212 
END IF 

WRITE (13,1032) 

DO 220 I = 1, NPRINT 
XX = (I-l)/FLOAT(NPRINT-l) 

CALL PRESS(XX,U,CP) 

PRINT 1040, XX, CP 
WRITE (14,1040) XX, CP 
220 WRITE (13,1040) XX, CP 
CLOSE (UNIT = 13) 

CLOSE (UNIT=14) 

C CALL LIBRARY ROUTINE TO CLEAR THE SCREEN, THEN PRINT HEADER 
CALL CLRSCRN 
PRINT * 

PRINT *, ' PROGRAM DUBLET RESULTS HAVE BEEN WRITTEN TO FILES: ' 
PRINT * 

PRINT *, ' DUBLET.DAT : DOUBLET STRENGTH DISTRIBUTION' 

PRINT *, ' SHAPE.DAT : BODY SURFACE COORDINATES' 

PRINT*, ' PRESSURE.DAT: SURFACE PRESSURE DISTRIBUTION' 

PRINT * 



PRINT * 

PRINT *, 'WOULD YOU LIKE TO PRINT THE RESULTS (Y/N)?' 
PRINT * 

READ 1002, PRINT 

IF (PRINT. EQ. 'Y'. OR. PRINT. EQ. 'y')THEN 
PRINT * 

PRINT *, 'WHICH OF THE FOLLOWING FILES DO YOU WANT?' 
PRINT * 

PRINT*, ' 1) DUBLET.DAT' 

PRINT*, ' 2) PRESSURE.DAT' 

PRINT*, ' 3) SHAPE.DAT' 

PRINT *, ' OR 4) ALL THREE FILES' 

PRINT * 

PRINT *, 'INPUT OPTION NO. (1,2,3, OR 4)' 

12 READ 1006, PRINTOPT 

IF (PRINTOPT .LT. 1 .OR. PRINTOPT . GT. 4) THEN 
PRINT *, 'INVALID ENTRY, ENTER INTEGER BETWEEN' 
PRINT *, 'ONE(l) AND F0UR(4). ' 

PRINT *, ' ' 

GO TO 12 
END IF 
ENDIF 

IF (PRINTOPT .EQ. 1) THEN 

CALL LIB$SPAWN( 'PRINT DUBLET.DAT') 

ENDIF 

IF (PRINTOPT .EQ. 2) THEN 

CALL LIB$SPAWN( 'PRINT PRESSURE.DAT') 

ENDIF 

IF (PRINTOPT .EQ. 3) THEN 

CALL LIB$SPAWN( 'PRINT SHAPE.DAT') 

ENDIF 

IF (PRINTOPT .EQ. 4) THEN 
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CALL LIB$SPAWN( 'PRINT DUBLET. DAT, PRESSURE. DAT, SHAPE. DAT' ) 
END IF 
PRINT * 

PRINT * 

PRINT *, 'WOULD YOU LIKE TO GRAPH THE RESULTS (Y/N)?' 

PRINT * 

READ 1002, GRAPH 

IF (GRAPH. EQ. ' Y' , OR. GRAPH. EQ. 'y')THEN 



PRINT 


* 






PRINT 


* 

y 


'WHICH OF THE FOLLOWING DATA FILES 


PRINT 


* 

y 


'DO YOU WANT TO 


GRAPH? ' 


PRINT 


* 






PRINT 


* 

y 


' 1) 


DUBLET. DAT' 


PRINT 


* 


' 2) 


PRESSURE.DAT' 


PRINT 


* 

y 


' 3) 


SHAPE. DAT' 


PRINT 


* 

y 


’ 4) 


NONE' 



PRINT * 

PRINT *, 'INPUT OPTION NO. (1,2,3 OR 4)' 

616 READ 1006, GRAPHOPT 

IF (GRAPHOPT .LT. 1 .OR. GRAPHOPT . GT. 4) THEN 
PRINT *, 'INVALID ENTRY, ENTER INTEGER BETWEEN' 

PRINT *, 'ONE(l) AND F0UR(4). ' 

PRINT ' 

GO TO 616 
ENDIF 

IF (GRAPHOPT .EQ. 1) THEN 
CALL GRAPH 1(NTYPE,XMAXY,TAU) 

C GET A HARDCOPY OF THIS GRAPHIC 

CALL LIB$SPAWN( 'RENDER/DEVICE=LA210/DRAFT_QUALITY/PAPER 
+SIZE=A Pl.UIS') 

CALL LIB$SPAWN( 'CONTINUE' ) 

PRINT *, ' ' 

PRINT *, 'WOULD YOU LIKE A PRINT OF THIS PLOT? (Y/N)' 
PRINT ' 

READ 1002, PLOTl 

IF (PLOTl. EQ. ' Y'. OR. PLOTl. EQ. 'y')THEN 
CALL LIB$SPAWN( 'PRINT PI. REN' ) 

ENDIF 
GO TO 46 
ENDIF 

IF (GRAPHOPT .EQ. 2) THEN 
. CALL GRAPH2(NTYPE,XMAXY,NPRINT,TAU,N) 

CALL LIB$SPAWN( 'RENDER/DEVICE=LA210/DRAFT_QUALITY/PAPER 
+SIZE=A P2.UIS') 

PRINT *, ' ' 

CALL LIB$SPAWN( 'CONTINUE' ) 

PRINT *, 'WOULD YOU LIKE A PRINT OF THIS PLOT? (Y/N)' 
PRINT *, ' ' 

READ 1002, PL0T2 

IF (PL0T2.EQ. 'Y'. OR. PL0T2. EQ. 'y')THEN 
CALL LIB$SPAWN(' PRINT P2. REN') 

ENDIF 
GO TO 46 
ENDIF 

IF (GRAPHOPT .EQ. 3) THEN 
CALL GRAPH3(NTYPE,XMAXY,N,TAU,DATP0) 
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CALL LIB$SPAWN( 'RENDER/DEVICE=LA210/DRAFT_QUALITY/PAPER_ 

+SIZE=A P3.UIS') 

PRINT ' 

CALL LIB$SPAWN( 'CONTINUE') 

PRINT *, 'WOULD YOU LIKE A PRINT OF THIS PLOT? (Y/N)' 

PRINT ' 

READ 1002, PLOT3 

IF (PL0T3.EQ. 'Y'.OR. PLOT3.EQ. 'y')THEN 
CALL LIB$SPAWN( 'PRINT P3. REN'' ) 

ENDIF 
GO TO 46 
ENDIF 

IF (GRAPHOPT .EQ. 4) THEN 
GO TO 64 
ENDIF 
ENDIF 

C OPTION TO MAKE ANOTHER RUN 

64 PRINT * 

PRINT *, ' DO YOU WISH TO: ' 

PRINT ' 1) MAKE ANOTHER RUN OR' 

PRINT ' 2) END THIS SESSION' 

PRINT *, ' ENTER 1 OR 2. ' 

PRINT * 

CALL QUERY (NANS) 

CALL CLRSCRN 

IF (NANS .EQ. 1) GO TO 10 

999 STOP 

1000 FORMAT(Al) 

1002 FORMAT(Al) 

1006 FORMAT(Il) 

1010 FORMAT(/,' DOUBLET STRENGTH DISTRIBUTION',/, 

+ ' M = M(I) FOR T(I) . LT. T . LT. T(I + 1)',//, 

+ 5X,'T(I)' ,5X,'M(I)/2',/) 

1012 FORMAT(/,9X, ' DOUBLET STRENGTH DISTRIBUTION',//, 

+ 14X,'T(I)' ,5X,'M(I)/2' ,/) 

1020 FORMAT(//,9X, ' BODY SHAPE - UPPER SURFACE ',//, 15X ,' X' , 9X, ' Y' , /) 
1030 FORMAT(//,' BODY SURFACE PRESSURE DISTRIBUTION',//, 

+ 6X, 'X' ,8X, 'CP' ,//, ' INPUT NUMBER OF PRESSURE COEFFICIENT', 

+ ' OUTPUT POINTS') 

1032 FORMAT(//,9X, ' BODY SURFACE PRESSURE DISTRIBUTION',//, 

+ 16X, 'X' ,8X, 'CP' ,//) 

1040 F0RMAT(9X,2F10. 4) 

END 



SUBROUTINE CLRSCRN 

LIBRARY ROUTINE TO CLEAR THE SCREEN. 

ISTAT = LIB$ERASE_PAGE (1,1) 

RETURN 
END 
C 

SUBROUTINE QUERY(NANS) 

C 

C ROUTINE TO TRAP ERRORS CAUSED BY IMPROPER RESPONSES TO QUESTIONS. 
C THE COMPUTER GENERATES AND ERROR WHEN A CHARACTER IS SUPPLIED TO 
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A QUESTION EXPECTING AN INTEGER OR REAL VALUE. 

NQTEST=0 
1 CONTINUE 

IF (NQTEST .GT. 0) THEN 

PRINT *, ’ CHARACTER VALUES ARE NOT VALID. ' 

PRINT *, ’ PLEASE ENTER AN INTEGER VALUE. ' 

END IF 

NQTEST = NQTEST + 1 
READ (5,*,ERR=1)NANS 
RETURN 
END 

SUBROUTINE FINDM (T,M,N,XS,XF) 

FIND DOUBLET STRENGTH TO MEET 
FLOW TANGENCY CONDITION 

REAL*4 T(100),M(100),XS,XF 
INTEGER N,R 

COMMON /COF/ A(101,111),NEQNS 
PI = 3. 1415926585 

NP = N + 1 

DO 100 I = 1,NP 

COSINE SPACING SCHEME FROM XS TO XF 
FRACT = .5*(1. - C0S(PI*(I-1)/FL0AT(N))) 

100 T(I) = XS + (XF - XS)*FRACT 

SET UP LINEAR SYSTEM OF EQUATIONS 

DO 210 I = 1,N 

XI = . 5*(T(I) + T(I+1)) 

YI = Y(XI) 

FACl = ATAN2(T(1) - XI, YI) 

DO 200 J = 1,N 

FAC2 = ATAN2(T(J+1) - XI, YI) 

A(I,J) = (FAC2 - FAC1)/YI 

FACl = FAC2 

200 CONTINUE 

A(I,NP) = 1.0 
210 CONTINUE 

SOLVE FOR DOUBLET STRENGTH 

NEQNS = N 

CALL GAUSS(l) 

DO 300 I = 1,N 

300 M(I) = A(I,NP) 

RETURN 

END 

SUBROUTINE FIX(VALMAX,VALMIN) 

C ARRAY = THE ARRAY WHICH IS BEING SORTED INTO ASCENDING ORDER 
C NUMBER= THE NUMBER OF ELEMENTS IN THE ARRAY TO BE SORTED 
C VALMAX= LARGEST VALUE IN THE ARRAY 
C VALMIN= SMALLEST VALUE IN THE ARRAY 
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REAL VALMAX.VALMIN 

INTEGER NUMBER 

LOGICAL SORTED 

COMMON /JACKEL/YPLOT,NP 

DIMENSION ARRAY(100),YPLOT(100) 

SORTED = .FALSE. 

NUMBER = NP 
DO 20 I = 1, NUMBER 
ARRAY(I) = YPLOT(I) 

20 CONTINUE 

30 IF (.NOT. SORTED) THEN 
SORTED = .TRUE. 

DO 40 I = 1, NUMBER - 1 

IF( ARRAY( I ) . GT. ARRAY( I+l) )THEN 
VALUE = ARRAY(I) 

ARRAY(I) = ARRAY(I+1) 

ARRAY(I+1) = VALUE 
SORTED = .FALSE. 

END IF 

40 CONTINUE 
GO TO 30 
END IF 

VALMAX = ARRAY( NUMBER) 

VALMIN = ARRAY(l) 

* THE FOLLOWING FILE IS CREATED TO CHECK THE VALIDITY OF THIS ROUTINE 
OPEN (UNIT=26,FILE='ARRAY3.DAT' , STATUS= ' NEW ’ ) 

DO 45 I = 1, NUMBER 

WRITE (17,55)ARRAY(I) 

45 CONTINUE 

WRITE ( 17, 65) VALMAX, VALMIN, NUMBER 
55 F0RMAT(1X,E11.4) 

65 FORMATC/, IX, 'VALMAX = ’ ,F6. 4,/, 'VALMIN = ', Ell. 4,/, 'NUMBER = ',13) 
CLOSE (UNIT=26) 

RETURN 

END 

SUBROUTINE GAUSS (NRHS) 

SOLUTION OF LINEAR ALGEBRAIC SYSTEM BY 
GAUSS ELIMINATION WITH PARTIAL PIVOTING 

A = COEFFICIENT MATRIX 

NEQNS = NUMBER OF EQUATIONS 

NRHS = NUMBER OF RIGHT HAND SIDES 

RIGHT-HAND SIDES AND SOLUTIONS STORED IN 
COLUMNS NEQNS+1 THRU NEQNS+NRHS OF ®A 

COMMON DX,DY,AR,PI 
COMMON /COF/ A(350,351) ,NEQNS 
NP = NEQNS + 1 

NTOT = NEQNS + NRHS 

GAUSS REDUCTION 
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DO 150 I = 2.NEQNS 



-- SEARCH FOR LARGEST ENTRY IN (I-l)TH COLUMN 
ON OR BELOW MAIN DIAGONAL 

IM =1-1 

IMAX = IM 

AMAX = ABS(A(IM,IM)) 

DO 110 J = I,NEQNS 

IF (AMAX .GE. ABS(A( J,IM))) GO TO 110 
IMAX = J 

AMAX = ABS(A(J,IM)) 
no CONTINUE 

SWITCH (I-l)TH AND IMAXTH EQUATIONS 

IF (IMAX .NE. IM) GO TO 140 
DO 130 J = IM,NT0T 
TEMP = A(IM,J) 

A(IM,J) = A(IMAX,J) 

A(IMAX,J) = temp 
130 CONTINUE 

ELIMINATE (I-l)TH UNKNOWN FROM 
ITH THRU (NEQNS)TH EQUATIONS 

140 DO 150 J = I,NEQNS 

R = A( J,IM)/A(IM,IM) 

DO 150 K = I,NTOT 

150 A(J,K) = A(J,K) - R*A(IM,K) 

BACK SUBSTITUTION 



DO 220 K = NP,NTOT 

A(NEQNS,K) = A(NEQNS,K)/A(NEQNS,NEQNS) 
DO 210 L = 2,NEQNS 

I = NEQNS + 1 - L 

IP =1+1 

DO 200 J = IP, NEQNS 

200 A(I,K) = A(I,K) - A(I,J)*A(J,K) 

210 A(I,K) = A(I,K)/A(I,I) 

220 CONTINUE 
RETURN 
END 



SUBROUTINE GRAPH1(NTYPE ,XMAXY,TAU) 

DEFINE I PACK ARRAY FOR LEGEND 
INTEGER*4 IPACK(35) 

REAL*4 T( 100) ,M( 100) ,XS ,XF,TAU,XMAXY,MIN,MAX 
INTEGER N,R,NTYPE,NP 
COMMON /MAIN/T,M,N,XS,XF 
COMMON /JACKEL/YPLOT,NP 
CHARACTER*40 LI 
DIMENSION YPLOT(IOO) 

C DEFINE AND ASSIGN LEGEND CHARACTER STRINGS 
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LI = 'DOUBLET STRENGTH^' 

C INITIALIZE THE GRAPHICS SYSTEM 
CALL INIT 

C LABEL X AND Y AXES USING SELF COUNTING STRINGS 
CALL XNAME( 'X' ,1) 

CALL YNAME( ' STRENGTH' ,8) 

C DEFINE PLOT AREA TO BE 6 INCHES BY 8 INCHES 
CALL AREA2D(6. 0,8. 0) 

C DEFINE HEADING LABEL 

CALL HEADINC ' DOUBLET STRENGTH DISTRIBUTION? ', -100,2. , 1) 

C PLOT ADDITIONAL TICK MARKS 

CALL XTICKS(l) 

CALL YTICKS(l) 

C PACK LEGEND LABELS INTO ARRAY IPACK 
CALL LINES(L1, IPACK, 1) 

C COSINE SPACING SCHEME FROM XS TO XF 
PI = 3. 1415926585 

NP = N + 1 

DO 100 I = 1,NP 

FRACT = ,5*(1. - C0S(PI*(I-1)/FL0AT(N))) 

100 T(I) = XS + (XF - XS)*FRACT 

C CREATE THE RESPECTIVE VALUES FOR YPLOT 

DO 207 I = 1,N+1 

YPLOT(I) = REAL(M(I)*3. 1415926585) 

207 CONTINUE 

CALL FIX(MAX,MIN) 

C SET UP AXIS 

CALL GRAF(0. ,.2,1. ,(MIN-. 1) ,. 05,(MAX+. 2)) 

C FRAME THE SUBPLOT AREA 

CALL FRAME 

C PLOT DUBLET STRENGTH DATA WITH A THICK LINE AND MARKER 15 
CALL MARKER(15) 

CALL THKCRV(.04) 

CALL CURVE (T, YPLOT, NP,1) 

C PLOT MESSAGES 

IF (NTYPE. EQ. 1) THEN 

CALL MESSAG( 'ELLIPTICAL AIRFOIL DOUBLET DISTRIBUTION? ', 100 , 
+ .5,6. ) 

CALL MESSAGC THICKNESS RATIO (TAU) = ? ' , 100 , . 5 ,5. 5) 

CALL REALN0(TAU,2,4. ,5.5) 

CALL MESSAGC 'NUMBER OF INTERVALS USED = ? ' , 100 , . 5 , 5. ) 

CALL INTNO(N, 'ABUT' ,' ABUT') 

END IF 

IF (NTYPE. EQ. 2) THEN 

CALL MESSAGC 'SYMMETRIC AIRFOIL DOUBLET DISTRIBUTION? ', 100 , 

+ .75,2.5) 

CALL MESSAGC 'THICKNESS RATIO (TAU) = ? ' , 100 , . 75 , 2. ) 

CALL REALN0(TAU,2,4. ,2. ) 

CALL MESSAGC 'MAXIMUM THICKNESS AT X = ? ' , 100 , . 75 , 1. 5) 

CALL REALN0(XMAXY,2,4. 1,1. 5) 

CALL MESSAGC 'NUMBER OF INTERVALS USED = ? ' , 100 , . 75 , 1. ) 

CALL INTNOCN, 'ABUT' , 'ABUT' ) 

ENDIF 

IF (NTYPE. EQ. 3) THEN 

CALL MESSAGC 'ARBITRARY SHAPE DOUBLET DISTRIBUTION?' 

+ ,100,. 75,1. 5) 
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CALL MESSAGC 'NUMBER OF INTERVALS USED = $ ' , 100 , . 75 , 1. ) 

CALL INTNO(N, 'ABUT' , 'ABUT' ) 

ENDIF 

C PLOT LEGEND 

CALL MYLEGNC 'DOUBLET STRENGTH? ', 100) 

C PLOT LEGEND 

CALL LEGEND(IPACK,1,3. 0,7. 0) 

C END PLOT 

CALL ENDPL(O) 

C CREATE GRAPHICS METAFILE PI. UIS 

CALL METAFL(l) 

C TERMINATE PLOT AT END OF PLOTTING SESSION 

CALL DONEPL 
RETURN 
END 

SUBROUTINE GRAPH2( NTYPE , XMAXY , NPRINT , TAU , N) 

DEFINE I PACK ARRAY FOR LEGEND 
INTEGER*4 IPACK(35) 

INTEGER NUM, NPRINT, NTYPE, N 
REAL XX( 100) ,CP( 100) , MAX, MIN, TAU, XMAXY 
CHARACTER*40 LI 
COMMON /ABLE/CP, NUM 
C READ ELEMENTS OF UNIT 14 INTO ARRAYS TO PLOT 
OPEN( UNIT=14 , FI LE= ' PRESS. DAT ' , STATUS= ' OLD ' ) 

DO 25 I = 1, NPRINT 
READ (14,*)XX(I) ,CP(I) 

25 CONTINUE 

NUM = NPRINT 
CL0SE(UNIT=14) 

CALL SCALER2(MAX,MIN) 

C DEFINE AND ASSIGN LEGEND CHARACTER STRINGS 

LI = 'CP DISTRIBUTION?' 

C INITIALIZE THE GRAPHICS SYSTEM 

CALL INIT 

C LABEL X AND Y AXES USING SELF COUNTING STRINGS 
CALL XNAME( 'X?' ,100) 

CALL YNAME( 'CP?' ,100) 

C DEFINE PLOT AREA TO BE 6 INCHES BY 8 INCHES 
CALL AREA2D(6. 0,8. 0) 

C DEFINE HEADING LABEL 

CALL HEADINCCP DISTRIBUTION? ', -100 ,2. ,1) 

C PLOT ADDITIONAL TICK MARKS 

CALL XTICKS(l) 

CALL YTICKS(l) 

C PACK LEGEND LABELS INTO ARRAY I PACK 
CALL LINES ( LI, IPACK,1) 

C SET UP AXIS 

CALL GRAF(0. 0,0. 2,1. 0,(MIN-. 1) ,((MAX-MIN)/5. ),(MAX+. 1)) 

C FRAME THE SUBPLOT AREA 

CALL FRAME 

C PLOT PRESSURE DISTRIBUTION DATA WITH A THICK LINE AND MARKER 15 
CALL MARKER(15) 

CALL THKCRV(.04) 
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CALL CURVE(XX,CP,NUM,1) 

C PLOT MESSAGES 

IF (NTYPE. EQ. 1) THEN 

CALL MESSAGC 'ELLIPTICAL AIRFOIL CP DISTRIBUTION? 100 , 

+ . 75,4. ) 

CALL MESSAGC 'THICKNESS RATIO (TAU) = $ ' , 100 , . 75 , 3. 5) 
CALL REALNO(TAU,2,4. ,3.5) 

CALL MESSAGC 'NUMBER OF INTERVALS USED = $ ' , 100 , . 75 , 3. 0) 
CALL INTNOCN, 'ABUT' ,' ABUT') 

END IF 

IF CNTYPE.EQ. 2) THEN 

CALL MESSAGC 'SYMMETRIC AIRFOIL CP DISTRIBUTION? ', 100 , 

+ . 75,6. 0) 

CALL MESSAGC 'THICKNESS RATIO CTAU) = ? ' , 100 , . 75 ,5. 5) 
CALL REALNOCTAU,2,4. 1,5. 5) 

CALL MESSAGC 'MAXIMUM THICKNESS AT X = ? ' , 100 , . 75 ,5. ) 
CALL REALN0CXMAXY,2,4. 1,5. ) 

CALL MESSAGC 'NUMBER OF INTERVALS USED = ? ' , 100 , . 75 ,4. 5) 
CALL INTNOCN, 'ABUT' , 'ABUT' ) 

ENDIF 

IF CNTYPE.EQ. 3) THEN 

CALL MESSAGC 'ARBITRARY SHAPE CP DISTRIBUTION?' 

+ ,100,. 75,5. 5) 

CALL MESSAGC 'NUMBER OF INTERVALS USED = ? ' , 100 , . 75 ,5. 0) 
CALL INTNOCN, 'ABUT' , 'ABUT' ) 

ENDIF 

C CHANGE LEGEND NAME TO "CP DISTRIBUTION" 

CALL MYLEGNC'CP DISTRIBUTION? ', 100) 

C PLOT LEGEND 

CALL LEGENDCIPACK,1,2. 0,7. 0) 

C END PLOT 

CALL ENDPLCO) 

C CREATE GRAPHICS METAFILE P2.UIS 

CALL METAFLC2) 

C TERMINATE PLOT AT END OF PLOTTING SESSION 

CALL DONEPL 
RETURN 
END 

SUBROUTINE GRAPH3 C NTYPE , XMAXY , N , TAU , DATPO ) 

DEFINE IPACK ARRAY FOR LEGEND 
INTEGER*4 IPACKC35) 

INTEGER NUM, NTYPE, N, DATPO 
REAL XXC100),YYC100),MAX,MIN,TAU,XMAY 
CHARACTER*40 LI 
COMMON / JACK/ YY, NUM 
C READ ELEMENTS OF UNIT 15 INTO ARRAYS TO PLOT 

OPENCUNIT=15,FILE='SHAPEBODY. DAT' , STATUS= ' OLD ' ) 

IF C NTYPE .LE. 2) THEN 
XXCl) = 0.0 
YYCl) = 0.0 
DO 25 I = 2,N+2 

READC15,*)XXCI),YYCI) 

25 CONTINUE 

XXCN+3) = 1. 0 
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YY(N+3) = 0. 0 
NUM = N + 3 
END IF 

IF (NTYPE .EQ. 3) THEN 
DO 35 I = l,DATPO+2 

READ(15,*)XX(I),YY(I) 

35 CONTINUE 

NUM = DATPO +2 
END IF 

CLOSE(UNIT=15) 

C CALL SCALER TO FIND THE MAX AND MIN VALUES OF THE Y ORDINATE ARRAY 
CALL SCALER(MAX,MIN) 

C DEFINE AND ASSIGN LEGEND CHARACTER STRINGS 
LI = 'AIRFOIL SHAPE$' 

C INITIALIZE THE GRAPHICS SYSTEM 

CALL INIT 

C LABEL X AND Y AXES USING SELF COUNTING STRINGS 
CALL XNAMEC ’X$' ,100) 

CALL YNAME( 'Y$' ,100) 

C DEFINE PLOT AREA TO BE 6 INCHES BY 8 INCHES 
CALL AREA2D(6. 0,8. 0) 

C DEFINE HEADING LABEL 

CALL HEADINC 'AIRFOIL SHAPE$ ' , - 100 , 2. , 1) 

C PLOT ADDITIONAL TICK MARKS 

CALL XTICKS(l) 

CALL YTICKS(l) 

C PACK LEGEND LABELS INTO ARRAY I PACK 

CALL LINES(L1,IPACK,1) 

C SET UP AXIS 

CALL GRAF(0. ,.2,1. ,0. ,( (MAX-MIN+. 4)/5. ),(MAX+. 4)) 

C FRAME THE SUBPLOT AREA 

CALL FRAME 

C PLOT PRESSURE DISTRIBUTION DATA WITH A THICK LINE AND MARKER 15 
CALL MARKER(15) 

CALL THKCRV(.04) 

CALL CURVE ( XX, YY, NUM, 1) 

C PLOT MESSAGES 

IF (NTYPE. EQ. 1) THEN 

CALL MESSAGC 'ELLIPTICAL AIRFOIL SHAPE? ',100, 

+ 1. ,5. ) 

CALL MESSAGC 'THICKNESS RATIO (TAU) = $ ' , 100, 1. ,4. 5) 

CALL REALN0(TAU,2,4. 5,4. 5) 

CALL MESSAGC 'NUMBER OF INTERVALS USED = $',100,1. ,4.0) 

CALL INTNOCN,' ABUT', 'ABUT') 

END IF 

IF (NTYPE. EQ. 2) THEN 

CALL MESSAGC 'SYMMETRIC AIRFOIL SHAPE? ',100, 

+ 1. ,5.0) 

CALL MESSAGC 'THICKNESS RATIO (TAU) = ?' , 100, 1. ,4. 5) 

CALL REALN0(TAU,2,4. 3,4. 5) 

CALL MESSAGC 'MAXIMUM THICKNESS AT X = ? ' , 100 , 1. ,4. ) 

CALL REALN0(XMAXY,2,4. 4,4. ) 

CALL MESSAGC 'NUMBER OF INTERVALS USED = ? ' , 100 , 1. ,3. 5 ) 

CALL INTNOC N ,' ABUT ',' ABUT ' ) 

ENDIF 

IF (NTYPE. EQ. 3) THEN 
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CALL MESSAGC 'ARBITRARY SHAPE$ ' , 100 , 1. ,4.5) 

CALL MESSAGC ' NUMBER OF INTERVALS USED = $',100,1. ,4.) 

CALL INTNOCN, 'ABUT' , 'ABUT' ) 

ENDIF 

C CHANGE LEGEND NAME TO "UPPER SURFACE ONLY" 

CALL MYLEGNC ' UPPER SURFACE$ ' , 100) 

C PLOT LEGEND 

CALL LEGEND(IPACK,1,3. 0,7. 0) 

C END PLOT 

CALL ENDPL(O) 

C CREATE GRAPHICS METAFILE P3. UIS 

CALL METAFLC3) 

C TERMINATE PLOT AT END OF PLOTTING SESSION 

CALL DONEPL 
RETURN 
END 

SUBROUTINE INTHV (NXTOL,NFTOL,NTYPE ,MAXIT, IOUT,UO ,U1) 

COMMON /MAIN/T,M,N,XS,XF 
DIMENSION T(100),M(100) 

SUBROUTINE TO FIND THE ROOTS OF f(x) = 0 USING THE 
INTERVAL HALVING METHOD 

IN THE PARAMETER LIST THE USER MUST PROVIDE: 

NXTOL = EXPONENT FOR X TOLERANCE VALUE 
NFTOL = EXPONENT FOR FUNCTION TOLERANCE VALUE 
NTYPE = SHAPE TYPE; ELLIPTICAL OR AIRFOIL 
MAXIT = MAXIMUM NUMBER OF ITERATIONS 
lOUT = 0 TO SUPPRESS ALL OUTPUT (TO DEVICE IW) 

1 TO OUTPUT FINAL RESULTS ONLY 

2 TO OUTPUT DETAILS FOR EACH ITERATION 
THE SUBROUTINE CALCULATES: 

XPREV, X = TWO INITIAL GUESSES, GIVEN N 
THE SUBROUTINE RETURNS: 

XS, XF = CURRENT X VALUES WHEN TERMINATION OCCURRED 
UO, U1 = CURRENT VELOCITY VALUES WHEN TERMINATION OCCURRED 
lEXIT = 1, 2, 3, 4 OR 7 (SEE FORMAT STATEMENTS 1 - 4 & 7) 

IW = 5 

XTOL = 10. **( -NXTOL) 

FTOL = 10. **( -NFTOL) 

CALCULATE INITIAL GUESS FOR XS AND XF, GIVEN N 
XS = 1. / FLOAT(N + 1) 

XSPREV = 10. **(-6) 

XF =1. - XS 
XFPREV =1. - XSPREV 

SET X VALUES FOR LEADING AND TRAILING EDGES FOR SUBROUTINE PRESS 
XLE = 0.0 
XTE = 1. 0 



ITERATE TO DETERMINE THE PROPER LOCATION FOR XF 

FIRST CHECK TO SEE THAT F(XF) & F(XFPREV) DIFFER IN SIGN 
SO THAT THE METHOD WILL CONVERGE. 
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EVALUATE PREVIOUS X VALUE 
CALL FINDM (T,M,N,XS,XFPREV) 

CALL PRESS (XTE,U1,CP) 

YFPREV = U1 

EVALUATE INITAL GUESS FOR X VALUE 
CALL FINDM (T,M,N ,XS ,XF) 

CALL PRESS (XTE,U1,CP) 

YF = U1 

IF (lOUT .GT. 1) WRITE (IW,5) XFPREV, YPREV, XF, YF 
IF (YFPREV*YF . GT. 0.0) THEN 
I = -2 
PRINT 201 
RETURN 
END IF 

COMPUTE SEQUENCE OF POINTS CONVERGING TO THE ROOT 
lEXIT = 1 
DO 10 K=l, MAXIT 

XR = (XFPREV + XF)/2. 0 

FOR THE ELLIPTIC CASE XS AND XF WILL BE EQUIDISTANT FROM THE EDGES. 
IF (NTYPE .LT. 2) THEN 
XS = ABS (1. - XR) 

END IF 

CALL FINDM (T,M,N,XS ,XR) 

CALL PRESS (XTE,U1,CP) 

Y = U1 

CHECK ON STOPPING CRITERIA 



DELTAXF = XFPREV-XR 
XERR = ABS(XFPREV-XR)/2. 0 

IF (lOUT .GT. 1) WRITE (IW,6) K,XR,Y,DELTA.XF 
IF (Y . EQ. 0. 0) lEXIT = 2 
IF (ABS(Y) .LE. FTOL) lEXIT = 3 
IF ( XERR . LE. XTOL) lEXIT = 7 
IF (lEXIT .GT. 1) GO TO 20 
IF (Y*YFPREV .GT. 0.0) THEN 
XFPREV = XR 
YFPREV = Y 
ELSE 

XF = XR 



YF = Y 
END IF 
10 CONTINUE 

THE MAXIMUM ITERATIONS HAS BEEN EXCEEDED, WITHOUT FINDING A ROOT. 
lEXIT = 4 



20 



EQ. 0) GO TO 30 



IF (lOUT 

IF (lEXIT .EQ. 1 
IF (lEXIT .EQ. 2 
IF (lEXIT .EQ. 3 
IF (lEXIT .EQ. 4 
30 CONTINUE 

FOR THE ELLIPTIC CASE XS ANND XF ARE DETERMINED, SO GO BACK. 



) WRITE (IW, 1) XR 

) WRITE (IW, 2) XR 

) WRITE (IW, 3) XR, NUMSIG 

) WRITE (IW, 4) MAXIT 



IF (NTYPE .LT. 2) THEN 

CALL FINDM (T,M,N,XS,XF) 
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CALL PRESS (XLE,U0,CP) 

GO TO 90 
END IF 

C NOW DO THE SAME FOR XS 

PRINT *, ' VALUE OBTAINED FOR XF' ,XF 
PRINT ' -- WORKING ON XS. ' 

C EVALUATE PREVIOUS X VALUE 

CALL FINDM (T,M,N,XSPREV,XF) 

CALL PRESS (XLE,U0,CP) 

YSPREV = UO 

C EVALUATE INITAL GUESS FOR X VALUE 
CALL FINDM (T,M,N,XS,XF) 

CALL PRESS (XLE,UO,CP) 

YS = UO 

IF (lOUT .GT. 1) WRITE (IW,5) XSPREV, YSPREV, XS , YS 
IF (YSPREV*YS .GT. 0.0) THEN 
I = -2 
PRINT 201 
RETURN 
END IF 

COMPUTE SEQUENCE OF POINTS CONVERGING TO THE ROOT 
lEXIT = 1 
DO 40 K=l, MAXIT 

XR = (XSPREV + XS)/2. 0 
CALL FINDM (T,M,N,XR,XF) 

CALL PRESS (XLE,UO,CP) 

Y = UO 

CHECK ON STOPPING CRITERIA 

DELTAXS = XSPREV-XR 
XERR = ABS(XSPREV-XR)/2. 0 

IF (lOUT .GT. 1) WRITE (IW,6) K,XR,Y, DELTAXS 
IF (Y . EQ. 0. 0) lEXIT = 2 
IF (ABS(Y) .LE. FTOL) lEXIT = 3 
IF ( XERR . LE. XTOL) lEXIT = 7 
IF (lEXIT .GT. 1) GO TO 50 
IF (Y*YSPREV .GT. 0.0) THEN 
XSPREV = XR 
YSPREV = Y 
ELSE 

XS = XR 
YS = Y 
END IF 



40 CONTINUE 



THE MAXIMUM 


ITERATIONS HAS 


BEEN 


EXCEEDED, WITHOUT FINDING A ROOT. 


lEXIT = 4 
















50 IF (lOUT . 


EQ. 


0) 


RETURN 








IF (lEXIT 


.EQ. 


1 


) 


WRITE 


(IW, 


1) 


XR 


IF (lEXIT 


.EQ. 


2 


) 


WRITE 


(iw, 


2) 


XR 


IF (lEXIT 


.EQ. 


3 


) 


WRITE 


(IW, 


3) 


XR, NUMSIG 


IF (lEXIT 


.EQ. 


4 


) 


WRITE 


(IW, 


4) 


MAXIT 


IF (lEXIT 


.EQ. 


7 


) 


WRITE 


(IW, 


7) 


XR, XTOL 



90 RETURN 
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THIS SHOULD RETURN WITH UO NEAR ZERO AND A GOOD VALUE OF XS. 

1 FORMATC ' OSLOPE = 0 WHEN X =',G12. 7,'. ITERATION DISCONTINUED. ') 

2 FORMATC ' OCOMPUTED F( ',G12. 7,') IS 0. ITERATION DISCONTINUED.') 

3 FORMATC ' OROOT: ' ,G12. 7 , ' . APPEARS TO BE ACCURATE TO ' , 1 1 , ' S. ' ) 

4 FORMATC 'ODES I RED ACCURACY IS NOT EVIDENT IN ',13,' ITERATIONS.') 

5 FORMATC ' OHALVING METHOD: XC-1! , XCO! ARE INITIAL GUESSES.',/, 

& '0 k' ,4X,'X = Xk' ,7X,'Y = FCX)' , 7X, ' X-XPREV' ,/, 

& ' -1 ', G12.7, E12.5, /,' 0 ', G12.7, E12.5) 

6 FORMATCI3, 3X, G12. 7, E12. 5, E15.5) 

7 FORMATC'OX LOCATION: ',G12. 7,' IS WITHIN X TOLERANCE OF ’,E12.5) 
201 FORMATC ' OFUNCTION HAS THE SAME SIGN AT BOTH INITIAL POSITIONS. ' 

& ,/,'OTHE BUILT-IN ITERATION SCHEME WILL NOT WORK, THEREFORE' 

& ,/,'OYOU MUST SELECT THE ENDPOINTS MANUALLY. ') 

END 



SUBROUTINE PRESS CX,U, CP) 

FIND PRESSURE COEFFICIENT CP AT CX,YCX)) 



100 



COMMON /MAIN/T,M,N,XS,XF 
DIMENSION TC100),MC100) 

REAL M 

YB = YCX) 

U =1.0 

V =0.0 

VFl = 1. /CCTCl) - X)**2 + YB*YB) 
UFl = CTCD - X)*VF1 

DO 100 J = 1,N 

VF2 = 1. /CCTCJ+1) - X)**2 + YB*YB) 
UF2 = CTCJ+1) - X)*VF2 

U = U + MCJ)*CUF2 - UFl) 

V = V - MCJ)*YB*CVF2 - VFl) 

VFl = VF2 

UFl = UF2 

CP =1.0- U*U - V*V 

RETURN 

END 



FUNCTION YCX) 



COMMON /FCN/ AX,TAU,NTYPE 
COMMON /DATA/COORXC 101) ,COORYC 101) 
COMMON /PROB/DATPO 
DIMENSION FDPCIOI) 

REAL FIRST, LATER, UP, DOWN, ARC 
INTEGER N,DATPO 



ORDINATE OF BODY CONTOUR 



IF CNTYPE .EQ. 1) THEN 

PROVIDE BODY ORDINATES FOR AN ELLIPSE OF THICKNESS RATIO TAU 
CCHORD HAS BEEN NONDIMENSIONALIZED, C=1.0) 

TO REDUCE THE NUMBER OF VARIABLES PASSED IN THE FUNCITON 
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STATEMENT, THE DUMMY VARIABLE AX PASSES TAU FOR THE ELLIPSOID 
CASE AND THE COEFFICIENT AX(TAU.XMAXY) FOR THE SYMMETRICAL 
AIRFOIL-LIKE CASE. 

Y = TAU * SQRT(X*(1. -X)) 

ELSEIF (NTYPE . EQ. 2) THEN 

PROVIDE BODY ORDINATES FOR A SYMMETRIC AIRFOIL-LIKE SHAPE 
(CHORD HAS BEEN NONDIMENSIONALIZED, C=1.0) 

Y = AX * SQRT(X)*(1-X) 

ELSE 

PROVIDE BODY ORDINATES FOR AN ARBITRARY BODY. TO DETERMINE 
THESE POINTS A CUBIC SPLINE INTERPOLATION SUBROUTINE WAS ADDED 
TO PROGRAM DUB LET. 

N = DATPO + 2 
XX = X 

CALL SPLINE(N,COORX,COORY,FDP) 

CALL SPEVALC N , COORX , COORY , FDP , XX , F ) 

Y = F 
END IF 
RETURN 
END 

SUBROUTINE SCALER( VALMAX, VALMIN) 

C ARRAY = THE ARRAY WHICH IS BEING SORTED INTO ASCENDING ORDER 
C NUMBER= THE NUMBER OF ELEMENTS IN THE ARRAY TO BE SORTED 
C VALMAX= LARGEST VALUE IN THE ARRAY 
C VALMIN= SMALLEST VALUE IN THE ARRAY 
REAL VALMAX, VALMIN 
INTEGER NUMBER 
LOGICAL SORTED 
COMMON /JACK/YY,NUM 
DIMENSION ARRAY(100),YY(100) 

SORTED = .FALSE. 

NUMBER = NUM 
DO 20 I = 1, NUMBER 
ARRAY(I) = YY(I) 

20 CONTINUE 

30 IF (. NOT. SORTED) THEN 
SORTED = .TRUE. 

DO 40 I = 1, NUMBER - 1 

IF(ARRAYd). GT. ARRAY( I+l) )THEN 
VALUE = ARRAY(I) 

ARRAY(I) = ARRAY(I+1) 

ARRAY(I+1) = VALUE 
SORTED = .FALSE. 

ENDIF 

40 CONTINUE 
GO TO 30 
ENDIF 

VALMAX = ARRAY( NUMBER) 
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VALMIN = ARRAY(l) 

* THE FOLLOWING FILE IS CREATED TO CHECK THE VALIDITY OF THIS ROUTINE 

OPEN (UNIT=16,FILE=’ ARRAY. DAT' ,STATUS=' NEW' ) 

DO 45 I = 1, NUMBER 

WRITE (16,55)ARRAY(I) 

45 CONTINUE 

WRITE (16, 65 )VALMAX, VALMIN, NUMBER 
55 F0RMAT(1X,E11. 4) 

65 F0RMAT(/,1X,'VALMAX = ' ,F6. 4,/, 'VALMIN = ', Ell. 4,/ ,' NUMBER = ',13) 
CLOSE (UNIT=16) 

RETURN 

END 

SUBROUTINE SCALER2(VALMAX, VALMIN) 

C 

C ARRAY = THE ARRAY WHICH IS BEING SORTED INTO ASCENDING ORDER 
C NUMBER= THE NUMBER OF ELEMENTS IN THE ARRAY TO BE SORTED 
C VALMAX= LARGEST VALUE IN THE ARRAY 
C VALMIN= SMALLEST VALUE IN THE ARRAY 
REAL VALMAX, VALMIN 
INTEGER NUMBER 
LOGICAL SORTED 
COMMON /ABLE/CP, NUM 
DIMENSION ARRAY(100),CP(100) 

SORTED = . FALSE. 

NUMBER = NUM 
DO 20 I = 1, NUMBER 
ARRAY(I) = CP(I) 

20 CONTINUE 

30 IF ( . NOT. SORTED) THEN 
SORTED = .TRUE. 

DO 40 I = 1, NUMBER - 1 

IF(ARRAYd). GT. ARRAY( I+l) )THEN 
VALUE = ARRAY(I) 

ARRAY(I) = ARRAY(I+1) 

ARRAY(I+1) = VALUE 
SORTED = .FALSE. 

ENDIF 

40 CONTINUE 
GO TO 30 
ENDIF 

VALMAX = ARRAY( NUMBER) 

VALMIN = ARRAY(l) 

* THE FOLLOWING FILE IS CREATED TO CHECK THE VALIDITY OF THIS ROUTINE 

OPEN (UNIT=17,FILE=' ARRAY2.DAT' , STATUS=' NEW' ) 

DO 45 I = 1, NUMBER 

WRITE (17,55)ARRAY(I) 

45 CONTINUE 

WRITE ( 17, 65)VALMAX, VALMIN, NUMBER 
55 F0RMAT(1X,E11.4) 

65 FORMATC/, IX, 'VALMAX = ' ,F6. 4,/, 'VALMIN = ', Ell. 4,/, 'NUMBER = ',13) 
CLOSE (UNIT=17) 

RETURN 

END 
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SUBROUTINE SPEVAL( N , COORX , COORY , FDP , XX , F ) 

THIS SUBROUTINE EVALUATES THE CUBIC SPLINE GIVEN 
THE DERIVATIVES COMPUTED BY SUBROUTINE SPLINE. 
THE INPUT PARAMETERS N,X,Y,FDP HAVE THE SAME 
MEANING AS IN SPLINE. 

XX = VALUE OF INDEPENDENT VARIABLE FOR WHICH 
AN INTERPOLATED VALUE IS REQUESTED 
F = THE INTERPOLATED RESULT 
DIMENSION COORXC 101) ,COORY( 101) ,FDP( 101) 

THE FIRST STEP IS TO FIND THE PROPER INTERVAL 
NMl = N -1 
DO 1 1=1, NMl 

IF (XX. LE. COORXC I+l)) GO TO 10 
1 CONTINUE 

NOW EVALUATE THE CUBIC 
10 DXM = XX - COORXC I) 

DXP = COORXC I+l) - XX 
DEL = COORXC I+l) - COORXC I) 

F = FDP(I)*DXP*(DXP*DXP/DEL - DEL)/6. 

1 +FDP(I+1)*DXM*(DXM*DXM/DEL - DEL)/6. 

2 +COORY(I)*DXP/DEL + COORY( I+1)*DXM/DEL 
RETURN 

END 



SUBROUTINE SPLINE (N, COORX, COORY, FDP) 

THIS SUBROUTINE COMPUTES THE SECOND DERIVATIVES NEEDED 
IN CUBIC SPLINE INTERPOLATION. THE INPUT DATA ARE: 

N = NUMBER OF DATA POINTS 

COORX = ARRAY CONTAINING THE VALUES OF THE INDEPENDENT VARIABLE 
(ASSUMED TO BE ASCENDING ORDER) 

COORY = ARRAY CONTAINING THE VALUES OF THE FUNCTION AT THE 
DATA POINTS GIVEN IN THE COORX ARRAY 
DIMENSION COORXC 101) ,COORY( 101) ,A( 101) ,B( 101) 

DIMENSION C(101),R(101),FDP(101) 

ALAMDA = 1 
NM2 = N - 2 
NMl = N - 1 

C(l) = COORXC2) - COORXC 1) 

DO 1 1=2, NMl 

C(I) = COORXC I+l) - COORXC I) 

A(I) = C(I-l) 

B(I) = 2.*(A(I) + CCD) 

R(I) = 6.*((COORY(I+l)-COORY(I))/C(I)-(COORY(I) 

+ -COORY(I-l))/C(I-l)) 

1 CONTINUE 

B(2) = B(2) + ALAMDA * C(l) 

B(NMl) = B(NMl) + ALAMDA * C(NMl) 

DO 2 1=3, NMl 
T = A(I)/B(I-1) 

B(I) = B(I) - T * C(I-l) 

R(I) = R(I) - T * R(I-l) 

2 CONTINUE 

FDP(NMl) = R(NM1)/B(NM1) 

DO 3 1=2, NM2 
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NMI = N - I 

FDP(NMI) = (R(NMI) - C(NMI) * FDP(NMI+1))/B(NMI) 
3 CONTINUE 

FDP(l) = ALAMDA * FDP(2) 

FDP(N) = ALAMDA * FDP(NMl) 

DESIRED DERIVATIVES HAVE NOW BEEN DETERMINED 

RETURN TO MAIN PROGRAM 

RETURN 

END 



FUNCTION YREF(XNUM) 

COMMON /LEROY/NUMERAL 

COMMON /BRAVO/NUMPTS 

COMMON /CHARLIE/NO 

COMMON /FLAGGER/FIGURE 

DIMENSION FDP( 101) ,XX( 101) ,YY( 101) 

DIMENSION XPOINTC 101) ,YPOINT( 101) ,XPOIN( 101) ,YPOIN( 101) 

NO = NUMPTS 

READ IN THE CURRENT SHAPE OF THE AIRFOIL 

IF(FIGURE. EQ. 2)NUMERAL=NUMERAL-2 
OPEN(UNIT=15,FILE='BODYSHAPE. DAT' , STATUS= ' OLD ' ) 

XX(1) = 0. 0 
YY(1) = 0. 0 
DO 30 I = 2,NUMERAL+1 
READ (15,*) XX(I),YY(I) 

30 CONTINUE 

XX(NUMERAL+2) = 1. 

YY(NUMERAL+2) = 0. 

CLOSE(UNIT=15) 

PROVIDE BODY ORDINATES FOR AN ARBITRARY . TO DETERMINE 
THESE POINTS A CUBIC SPLINE INTERPOLATION SUBROUTINE WAS ADDED 
TO PROGRAM NEW_PANEL. 

THE AIRFOIL SHAPE IS BEING SPLIT INTO UPPER AND LOWER SURFACES AND 
THEN FORMATTED TO BE USED WITH THE SPLINE/SPEVAL ROUTINES. 

NOB = INT(NUMPTS/2)+l 
DO I=l,INT(NUMPTS/2)+l 
DUMMY=XX(I) 

DUM =YY(I) 

XPOINTC I )=DUMMY 
YPOINT(I)=DUM 
END DO 

DO I=INT(NUMPTS/2)+2, NUMPTS 
DUMM=XX( I ) 

DU =YY(I) 

XPOIN(I)=DUMM 
YPOIN(I)=DU 
END DO 

XPOIN(NUMPTS+l)=l. 

YPOIN(NUMPTS+1)=0. 

CALL SORTNUM(XPOINT,YPOINT,NOB) 
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CALL S0RTNUM(XP0IN,YP0IN,N0B-1) 

UPPER SURFACE Y COORDINATE DETERMINATION 

IF (XNUM. GT. 0. )THEN 
N = INT(NUMPTS/2)+l 
XPT = XNUM 

CALL SPLINE(N,XPOINT,YPOINT,FDP) 

CALL SPEVALC N , XPO INT , YPO INT , FDP , XPT , F ) 
YREF = F 
END IF 

LOWER SURFACE Y COORDINATE DETERMINATION 

IF (XNUM. LT. 0. ) then 
N = INTCNUMPTS/2) 

XPT = XNUM 

CALL SPLINE(N,XPOIN,YPOIN,FDP) 

CALL SPEVALC N , XPOIN , YPO IN , FDP , XPT , F) 
YREF = F 
END IF 
RETURN 
END 
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APPENDIX G. PROGRAM NEVV_PANEL COMPUTER CODE 

PROGRAM NEW_PANEL 
C 

C *** MODIFIED FOR USE ON THE MICROVAX/2000 BY J. A. CAMPBELL (JUL 88) 

C UPDATES MADE BY C. M. MACALLISTER JAN-NOV 89 (CMM) 

Q *********Vc**Vfyf********':V****Vr*Vc**VfV-***y-*******Vf ********************** 

C PROGRAM NEW.PANEL 

C 

C SMITH-HESS (DOUGLAS) PANEL METHOD 

C FOR SINGLE -ELEMENT LIFTING AIRFOIL 

C IN TWO-DIMENSIONAL INCOMPRESSIBLE FLOW 

C 

C ORIGINAL IBM MAINFRAME PROGRAM WAS ADAPTED FROM JACK MORAN'S BOOK 
C 'AN INTRODUCTION TO THEORETICAL AND COMPUTATIONAL AERODYNAMICS', 

C WILEY AND SONS, NEW YORK 1984. THE LISTING IS FOUND ON PAGE 118. 

C 

C PROGRAM FLEXIBILITY AND USER INTERFACE WAS REVISED FOR 
C PROFESSOR J. V. HEALEY BY JOHN CAMPBELL. APRIL 1988. 

C ADDITIONAL PROGRAM UPDATES TO INCLUDE PRINTING ROUTINES, 

C PROCESSING CORRECTIONS, GRAPHICAL ANALYSIS, AND VISCOUS 

C INTERACTION ADDAPTATION WERE MADE BY CRAIG MACALLISTER 

C IN JAN-NOV 1989. (CMM) 

C 

C THE VISCOUS INTERACTION ADAPTATION WAS REALIZED USING A 

C FORTRAN PROGRAM DEVELOPED BY DR. T. CEBECI AND DR. H. 

C B. KELLER. THE ORIGINAL VERSION OF THIS PROGRAM IS 

C CURRENTLY AVAILABLE FOR USE ON THE IBM MAINFRAME AT THE 

C NAVAL POSTGRADUATE SCHOOL, ACCOUNT 4632P. IN ORDER TO 

C USE DR. CEBECI 'S PROGRAM IT IS NECESSARY TO INPUT THE 

C POTENTIAL FLOW SOLUTION OVER A SECTION SHAPE. IN PARTICU- 

C LAR, THE CP DISTRIBUTION OR THE VELOCITY DISTRIBUTION OVER 

C AS SURFACE IS REQUIRED. SUCH INFORMATION IS OBTAINED QUITE 

C READILY THROUGH THE EXECUTION OF THE NEW_PANEL PROGRAM. IN 

C FACT, THE CP DISTRIBUTION CREATED BY NEW_PANEL IS INTER- 

C ACTIVELY ADAPTED, SORTED AND INPUTTED TO THE CEBECI PROGRAM 

C UPON THE USER'S REQUEST. 

C 

C THIS PROGRAM PROVIDES THE BODY COORDINATES AND THE SURFACE 
C PRESSURE DISTRIBUTION ABOUT A SINGLE ELEMENT LIFTING AIRFOIL 
C IN TWO-DIMENSIONAL FLOW. 

C 

C ESTIMATED VALUES FOR LIFT COEFFICIENT AND THE MOMENT COEFFICIENT 
C ABOUT THE LEADING EDGE AND QUARTER CHORD ARE DETERMINED FROM THE 
C PRESSURE COEFFICIENTS OF EACH PANEL. 

C 

C YOU MAY PROVIDE ACTUAL AIRFOIL SURFACE COORDINATE DATA VALUES OR 
C HAVE THE COMPUTER GENERATE AN APPROXIMATION FOR THE COORDINATES 

C OF A NACA XXXX OR 230XX AIRFOIL SECTION. 

C 

C IF YOU DESIRE TO ENTER THE SURFACE COORDINATE VALUES, SEVERAL 
C OPTIONS ARE AVAILABLE. YOU MAY ENTER THEM 1) FROM A DATA FILE, 

C 2) FROM THE KEYBOARD OR 3) USING DATA STATEMENTS ALREADY ENTERED 

C AT THE END OF THE MAIN PROGRAM LISTING. 

C 

C IF INPUTTING YOUR OWN DATA, REMEMBER TO START AT THE TRAILING EDGE 
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(X/C = 1.0), AND WORK TOWARDS THE LEADING EDGE, ENTERING THE LOWER 
SIDE FIRST, FOLLOWED BY THE UPPER SURFACE. DO NOT ENTER THE 
TRAILING EDGE TWICE. TRY TO ENTER A SUFFICIENT NUMBER OF POINTS 
NEAR THE NOSE FOR GOOD RESOLUTION. 

*** NOTE: TO SATISFY THE KUTTA CONDITION, X VALUES FOR POINTS 
2 AND N MUST BE THE SAME. THIS ENSURES THAT THE LAST 
PANELS, UPPER AND LOWER, ARE OF EQUAL SIZE. 

CD IS JUST AN INDICATOR OF NUMERICAL ACCURACY OF THIS 
PROGRAM. VALUE OF CD SHOULD BE NEAR ZERO. 

IF USING DATA STMTS OR AN INPUT FILE, REMEMBER THE NUMBER 
OF DATA VALUES AS YOU WILL BE ASKED FOR THIS BY THE PROGRAM. 

USE OF THE DATA STATEMENTS REQUIRES THAT PROGRAM BE 
MODIFIED IN ADVANCE BY MOVING THEM TO THE CORRECT LOCATION. 

INTEGER NANS, FLAG, FIGURE 

REAL L I FTA , MOMENTA , MENTA 

DIMENSION XX( 101) , YCORD( 101) , YCOORD( 101) ,DLS( 101) ,THT(' 101) 
DIMENSION ANGLE( 13) ,CLA( 13) ,CMA( 13) ,CMB( 13) ,YY( 101) ,XREF( 101) 
DIMENSION YDAT(lOl) 

*** NOTE: IF YOU CHANGE SIZE OF X AND Y, CHANGE N BELOW ALSO| *** 

CHARACTER*! PRINT , GRAPH , PLOTl , PLOT2 , PLOT3 , PL0T4 , PLOTS 
CHARACTER*! VISCOUS , PRINTER 
INTEGER PRINTOPT , GRAPHOPT , THICKOPT , VISOPT 
DATA X, Y / 101*0. ,101*0. / 

DATA ANGLE/-8. ,-6. ,-4. ,-2. ,0. ,2. ,4. ,6. ,8. ,10. ,12. ,14. ,16. / 

COMMON /EXTRA/LIFTA , MOMENTA , MENTA 
COMMON /DATA/XC 10 1 ) , Y( 10 1 ) 

COMMON /BOD/ NODTOT,COSTHE( 100) ,SINTHE( 100) ,NFLAG 
COMMON /PAR/ NACA,TAU,EPSMAX,PTMAX 
COMMON /FLAGGER/FIGURE 
COMMON /FINAL/FLAG, XREF,YCORD 
COMMON /COF/ A(101,111),KUTTA 
COMMON /BRAVO/NUMPTS 
COMMON /NUM/ PI,PI2INV 
COMMON /CPD/ CP(IOO) 

COMMON /PEN/ CLA,CMA, ANGLE, CMB 

C IF USING DATA STMTS FOR X AND Y VALUES, PLACE LINES HERE. 

C Mric following data is for the NASA LS(1)-0013 AIRFOIL *** 

DATA NUPPER, NLOWER /14,14/ 

DATA (X( I), 1=1, 28)/ 1. 0,. 90,. 80,. 70,. 60,. 50,. 40,. 30,. 20,. 10, 

1 0. 07535,0. 05,0. 0247,0. 01255,0. 0,0. 01301,0. 02505,0. 04993,0. 07498, 

2 0. 10,. 20, < 30,. 40,. 50,. 60,. 70,. 80,. 90/ 

C *** NOTE: VALUE FOR TRAILING EDGE IS SET TO 0. 000 VS ACT THICKNESS * 
DATA (Y(I),I=l,28)/0. 00000,-. 01165,-. 02654,-. 04196,-. 05459, 

1 -. 06209,-. 06453,-. 06316,-. 05755,-. 04543,-. 04070,-. 03462, 

1 -.02612,-. 01938,0. 0,. 01892,. 02583,. 03465,. 04075,. 04541, 
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2 . 05750,. 06307,. 06432,. 06203,. 05446,-04183,. 02638,-01172/ 

Q Vc*****:?r******VfVcyc*****VcVfVc*Vf**Vc'5V')V***VfVf***ycyc*Vf*VoVycycycVr*Vf*VcVf*VfVc’'}VVf***'}V***-jV 

PI = 3. 1415926585 

C *** MAKE SURE THAT N CORRESPONDS TO THE SIZE OF X AND Y DIMENSION ** 

59 N = 100 
FLAG = 1 

Q '}V'>V**VfVf*Vc')VVc*Vc*Vfyc'>V**********yc******VfVc****Vc***Vr**Vf*Vc***Vf***yc********Vf'>V* 

c 

C OPEN FILE FOR BODY SURFACE COORDINATE OUTPUT 

60 OPEN (UNIT=11, 

2 FILE= 'PB0DY.DAT', 

2 0RGANIZATI0N= 'SEQUENTIAL' , 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE', 

2 FORM= 'FORMATTED', 

2 STATUS= 'NEW') 

C OPEN FILE FOR PRESSURE COEFFICIENT OUTPUT 
OPEN (UNIT=12, 

2 FILE= 'PPRESS.DAT', 

2 0RGANIZATI0N= 'SEQUENTIAL' , 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE' , 

2 F0RM= 'FORMATTED', 

2 STATUS= 'NEW') 

OPEN ANOTHER FILE FOR BODY SURFACE PRESSURE DISTRIBUTION OUTPUT 
OPEN (UNIT=14, 

2 FILE= 'PRESSER.DAT', 

2 ORGANIZATION= 'SEQUENTIAL' , 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE', 

2 FORM= 'FORMATTED', 

2 STATUS= 'NEW') 

OPEN ANOTHER FILE FOR BODY SHAPE OUTPUT 
OPEN (UNIT=15, 

2 FILE= 'BODYSHAPE.DAT', 

2 ORGANIZATION= 'SEQUENTIAL', 

2 ACCESS= 'SEQUENTIAL' , 

2 RECORDTYPE= 'VARIABLE', 

2 FORM= 'FORMATTED', 

2 STATUS= 'NEW') 



0PEN(UNIT=46,FILE='ME4. DAT' , STATUS= ' NEW ' ) 
DO 1= 1,NUMPTS 
WRITE(46,999)X(I),Y(I) 

END DO 

CLOSE(UNIT=46) 

CALL INDATA(X,Y,N,NLOWER,NUPPER) 

IF (FLAG. EQ. 2)THEN 
NLOWER = NUMPTS/2 
N = NUMPTS 
NUPPER = NUMPTS/2 

0PEN(UNIT=43,FILE='ME1.DAT' , STATUS=' NEW' ) 
DO 1= 1, NUMPTS 
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WRITE(43,999)X(I),Y(I) 

END DO 

CL0SE(UNIT=43) 

ENDIF 

CALL SETUP(X,Y,N,NLOWER,NUPPER) 

CHECK THE INPUT OF THE AIRFOIL SHAPE DATA( OPTIONAL) 

OPEN (66,FILE='MAKE.DAT’ , STATUS= ' NEW ' ) 

DO I = 1,NUMPTS 

WRITE(66,999)X(I),Y(I) 

END DO 

999 F0RMAT(1X,F8. 4,2X,F8. 4) 

CLOSE (UNIT=66) 

IF (FLAG. EQ. 2)G0 TO 908 
100 PRINT 1000 

READ (5,*) ALPHA 

IF (ALPHA .GT. 90.) GO TO 200 

ADA = ATPHA 

COSALF = C0S(ALPHA*PI/180. ) 

SINALF = SIN(ALPHA*PI/180. ) 

CALL COFI SH( SINALF , COSALF , X , Y , N , NLOWER , NUPPER) 

CALL GAUSS(l) 

OPEN(UNIT=45',FILE='ME3. DAT' , STATUS=' NEW' ) 

DO 1= 1,NUMPTS 
WRITE(45,999)X(I) ,Y(I) 

END DO 

CLOSE(UNIT=45) 

CALL VELD I S ( S INALF , COSALF , X , Y , N , NLOWER , NUPPER , ALPHA ) 

CALL FANDM( SINALF , COSALF , X , Y , N , NLOWER , NUPPER) 

00 CONTINUE 

CL0SE(UNIT=11) 

CL0SE(UNIT=12) 

CL0SE(UNIT=14) 

CL0SE(UNIT=15) 

908 IF (FLAG. EQ. 2)THEN 

CALL COFISH( SINALF , COSALF , X , Y , N , NLOWER , NUPPER) 

CALL GAUSS(l) 

0PEN(UNIT=44,FILE='ME2. DAT' , STATUS= ' NEW ' ) 

DO 1= 1,NUMPTS 
WRITE(44,999)X(I),Y(I) 

END DO 

CLOSE (UNIT=44) 

CALL VELDIS ( SINALF , COSALF , X , Y , N , NLOWER , NUPPER , ALPHA) 
CALL FANDM( SINALF , COSALF , X , Y , N , NLOWER , NUPPER) 
CLOSE(UNIT=ll) 

CLOSE(UNIT=12) 

CL0SE(UNIT=14) 

CL0SE(UNIT=15) 

ENDIF 

C CALL LIBRARY ROUTINE TO CLEAR THE SCREEN, THEN PRINT HEADER 
CALL CLRSCRN 
PRINT * 
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12 



46 



65 



PRINT *, ' PROGRAM PANEL RESULTS HAVE BE-EN WRITTEN TO FILES: ' 
PRINT * 

PRINT *, ' PBODY.DAT ; BODY SURFACE COORDINATES' 

PRINT *, ' PPRESS.DAT : SURFACE PRESSURE DISTRIBUTION' 

PRINT * 

PRINT * 

PRINT *, 'would you LIKE TO PRINT THE RESULTS (Y/N)?' 

PRINT * 

READ 1002, PRINT 

IF ( PRINT. EQ. ' Y' . OR. PRINT. EQ. ' y ' )THEN 
PRINT * 

PRINT *, 'WHICH OF THE FOLLOWING FILES DO YOU WANT?' 

PRINT * 



PRINT 


* » 

> 




1) 


PBODY. DAT' 


PRINT 


* * 

y 




2) 


PPRESS. DAT 


PRINT 


rk ^ 

y 


OR 


3) 


BOTH FILES 



PRINT * 

PRINT *, 'INPUT OPTION NO. (1,2, OR 3)' 

READ 1006, PRINTOPT 

IF (PRINTOPT .LT. 1 .OR. PRINTOPT . GT. 3) THEN 
PRINT *, 'INVALID ENTRY, ENTER INTEGER BETWEEN' 
PRINT *, 'ONE(l) AND THREE(3). ' 

PRINT *, ' ' 

GO TO 12 
END IF 
END IF 

IF (PRINTOPT .EQ. 1) THEN 

CALL LIB$SPAWN( 'PRINT PBODY.DAT') 

ENDIF 

IF (PRINTOPT .EQ. 2) THEN 

CALL LIB$SPAWN( 'PRINT PPRESS.DAT') 

ENDIF 

IF (PRINTOPT .EQ. 3) THEN 

CALL LIB$SPAWN( ' PRINT PBODY. DAT,PPRESS. DAT' ) 
ENDIF 

CALL CLRSCRN 
PRINT * 



PRINT * 

PRINT *, 'WOULD YOU LIKE TO GRAPH THE RESULTS (Y/N)?' 
PRINT * 

READ 1002, GRAPH 

IF (GRAPH. EQ. ' Y' . OR. GRAPH. EQ. ' y' )THEN 
CALL CLRSCRN 
PRINT * 

PRINT *, 'WHICH OF THE FOLLOWING DATA OUTPUTS' 

PRINT *, ' DO YOU WANT TO PLOT?' 

PRINT * 

PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 

PRINT *, 'INPUT OPTION NO. (1,2,3 OR 4)' 

READ 1006, GRAPHOPT 

IF (GRAPHOPT .LT. 1 .OR. GRAPHOPT . GT. 4) THEN 



* 

J 

* 

J 

* 

i 



1 ) 

2 ) 

3) 

4 ) 



PPRESS.DAT (CP DISTRIBUTION) ' 
PBODY.DAT (AIRFOIL SHAPE) ' 

CL VS. ANGLE OF ATTACK ' 

& CM C/4 VS. ANGLE OF ATTACK ' 
NONE' 



124 



PRINT *, 'INVALID ENTRY, ENTER INTEGER BETWEEN' 

PRINT *, 'ONE(l) AND F0UR(4). ' 

PRINT ' 

GO TO 65 
END IF 

IF (GRAPHOPT .EQ. 1) THEN 
CALL GRAFl 

C GET A HARDCOPY OF THIS GRAPHIC 

CALL LIB$SPAWN( 'RENDER/DEVICE=LA210/DRAFT_QUALITY/PAPER_ 
+SIZE=A Pl.UIS') 

PRINT ' 

PRINT *, 'WOULD YOU LIKE A PRINT OF THIS PLOT? (Y/N)' 
PRINT ' 

READ 1002, PLOTl 

IF (PLOTl. EQ. 'Y'. OR. PLOTl. EQ. 'y' )THEN 
CALL LIB$SPAWN( 'PRINT PI. REN') 

END IF 
GO TO 46 
ENDIF 

IF (GRAPHOPT .EQ. 2) THEN 
CALL GRAF 2 

CALL LIB$SPAWN( 'RENDER/DEVICE=LA210/DRAFT_ QUALITY/PAPER. 
+SIZE=A P2.UIS') 

PRINT ' 

PRINT *, 'WOULD YOU LIKE A PRINT OF THIS PLOT? (Y/N)' 
PRINT *, ' ' 

READ 1002, PLOT2 

IF (PLOT2. EQ. 'Y' . OR. PLOT2. EQ. 'y' )THEN 
CALL LIB$SPAWN( 'PRINT P2. REN') 

ENDIF 
GO TO 46 
ENDIF 

IF (GRAPHOPT .EQ. 3) THEN 

Q Vr Vr Vr Vc Vr Vr Vr Vr Vr V' ic Vr Vr Vc Vr 

OPEN (UNIT=11, 

2 FILE= 'PBODY.DAT', 

2 ORGANIZATION= 'SEQUENTIAL' , 

2 ACCESS= ' SEQUENTIAL' , 

2 RECORDTYPE= 'VARIABLE', 

2 FORM= 'FORMATTED', 

2 STATUS= 'UNKNOWN') 

OPEN (UNIT=12, 

2 FILE= 'PPRESS.DAT', 

2 ORGANIZATION= 'SEQUENTIAL' , 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE', 

2 FORM= 'FORMATTED', 

2 STATUS= 'UNKNOWN') 

OPEN (UNIT=14, 

2 FILE= 'PRESSER.DAT', 

2 ORGANIZATION= 'SEQUENTIAL' , 

2 ACCESS= 'SEQUENTIAL', 

2 RECORDTYPE= 'VARIABLE' , 

2 FORM= 'FORMATTED', 

2 STATUS= 'UNKNOWN') 

OPEN (UNIT=15, 
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C 



C 



68 



2 FILE= 'BODYSHAPE. DAT' , 

2 0RGANIZATI0N= 'SEQUENTIAL' , 

2 ACCESS= 'SEQUENTIAL' , 

2 RECORDTYPE= 'VARIABLE' , 

2 FORM= 'FORMATTED', 

2 STATUS= 'unknown') 

DO 45 I = 1,13 
ALPHA = ANGLE(I) 

COSALF = COS(ALPHA*PI/180. ) 

SINALF = SIN(ALPHA*PI/180. ) 

CALL COFISHC S INALF , COS ALF , X , Y , N , NLOWER , NUPPER) 

CALL GAUSS(l) 

CALL VELDIS ( S INALF , COSALF , X , Y , N , NLOWER , NUPPER , ALPHA) 
CALL FANDM( S INALF , COS ALF , X , Y , N , NLOWER , NUPPER) 

CLA(I) = LIFTA 
CMA(I) = MOMENTA 
CMB(I) = MENTA 
CLOSE(UNIT=ll) 

CL0SE(UNIT=12) 

CL0SE(UNIT=14) 

CL0SE(UNIT=15) 

CONTINUE 
END IF 

IF (GRAPHOPT .EQ. 3)THEN 
CALL GRAF3 

GET A HARDCOPY OF THIS GRAPHIC 

CALL LIB$SPAWN( 'RENDER/DEVICE=LA210/DRAFT_QUALITY/PAPER 
+SIZE=A P3.UIS') 

PRINT ' 

PRINT *, 'WOULD YOU LIKE A PRINT OF THIS PLOT? (Y/N)' 
PRINT *, ' ' 

READ 1002, PLOT3 

IF ( PL0T3. EQ. ' Y' . OR. PLOT3. EQ. ' y ' )THEN 
CALL LIB$SPAWN( ' PRINT P3. REN' ) 

ENDIF 
GO TO 46 
ENDIF 

IF (GRAPHOPT .EQ. 4 .AND. FLAG . EQ. 1) GO TO 68 
IF (GRAPHOPT .EQ. 4 .AND. FLAG . EQ. 2) GO TO 64 
ENDIF 

IF (FLAG. EQ. 2)G0 TO 64 



CALL CLRSCRN 
PRINT *, ' ' 

PRINT *, ' WOULD YOU LIKE TO ANALYZE VISCOUS EFFECTS' 
PRINT *, ' FOR THIS AIRFOIL (Y/N) ?' 

PRINT *, ' ' 

READ 1002, VISCOUS 

IF (VISCOUS. EQ. 'Y'. OR. VISCOUS. EQ. 'y')THEN 
FLAG = 2 



FIGURE = 2 
PRINT *, ' 
PRINT *, ' 
PRINT *, ' 
PRINT *, ' 
PRINT *, ' 



VISCOUS BOUNDRY LAYER ANALYSIS' 

I 

*** INPUT DATA OPTION ***' 

t 

WHAT INPUT SOURCE WOULD YOU LIKE TO USE?' 
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ooooooooooooooooooooooooo onoo 



1) DATA FILE "BL2D.DAT" OR ' 

2) NEW_PANEL CP DISTRIBUTION JUST CREATED' 

3) QUIT PROGRAM' 

f 

ENTER 1, 2, OR 3' 



PRINT *, 

PRINT *, 

PRINT *, 

PRINT *, 

PRINT *, 

PRINT *, 

PRINT *, 

42 READ *, VI SORT 

IF (VISOPT. LT. 1. OR. VISOPT. GT. 3)THEN 

PRINT ' INVALID ENTRY | TRY AGAIN| ' 

PRINT *, ' ENTER 1, 2, OR 3' 

GO TO 42 
END IF 

IF (VISOPT. EQ. 3) GO TO 1111 
CALL CIB(VISOPT) 

PRINT ■ 

THE BOUNDRY LAYER RESULTS HAVE BEEN 

WRITTEN TO FILE "BL2D. OUT'" 



WOULD YOU LIKE TO PRINT THESE RESULTS?' 

I 



PRINT *, 

PRINT *, 

PRINT *, 

PRINT 
PRINT 
READ 1002, PRINTER 

IF (PRINTER. EQ. ' Y' . OR. PRINTER. EQ. 'y')THEN 
CALL LIB $SPAWN( 'PRINT BL2D. OUT') 

ENDIF 

IF (VISOPT. EQ. 1)G0 TO 64 



THE FOLLOWING "GO TO" STATEMENT WAS PUT IN TO CIRCUMVENT 
PROCESSING THE PANEL METHOD AGAIN. 



GO TO 64 

NUMPTS = NUMPTS-2 

OPEN(UNIT=63,FILE='VISC. DAT’ , STATUS=’ UNKNOWN' ) 
READ(63,920)XP,DLSN,THTN 
READ(63,920)XP,DLSN,THTN 
READ(63,920)XP,DLSN,THTN 
DO I=l,INT(NUMPTS/2) 

READ(63,920)XREF(I) ,DLS( I) ,THT( I) 

XNUM = XREF(I) 

YYY = YREF(XNUM) 

YCORD(I)=YYY 
END DO 

DO I=(INT(NUMPTS/2)+l) , NUMPTS 

READ(63,920)XREF(I),DLS(I) ,THT(I) 

XNUM = XREF(I) 

YYY = YREF(XNUM) 

YCORD(I) = -YYY 
END DO 

CLOSE(UNIT=63) 

CHECKING THE DATA COORDINATES TO SEE IF IN PROPER ORDER 
OPEN (UNIT=77,FILE='RET.DAT' , STATUS= ' NEW ’ ) 

DO I=1,NU?JPTS 

WRITE (77,991)XREF(I) ,YCORD(I) 

'91 F0RMAT(1X,F8.4,2X,F8. 4) 

END DO 

CLOSE (UNIT=77) 
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69 



64 



PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 



WHICH TYPE OF BOUNDRY LAYER THICPCNESS ' 
WOULD YOU LIKE TO ANALYZE* 

I 

1) DISPLACEMENT THICKNESS* 

2) MOMENTUM THICKNESS* 

3) NONE--QUIT VISCOUS ANALYSIS* 



READ 1006 

IF (THICKOPT 

* ' 
> 

* ' 



THICKOPT 
LT. 1 

I 



* 

* 

69 



INVALID ENTRY, ENTER INTEGER BETWEEN* 
ONE(l) AND THREE(3). * 



OR. THICKOPT ,GT. 3) THEN 

PRINT 
PRINT 
PRINT 
PRINT 
GO TO 
ENDIF 

IF (THICKOPT. EQ. 1) THEN 
DO I=1,NUMPTS 
THICK= DLS(I) 

IF (YCORD(I). LT. 0. )YCOR = YCORD( I) -THICK 
IF (YCORD(I). GE. 0. )YCOR = YCORD( I )+THICK 
YCORD(I) = YCOR 
END DO 
ENDIF 

IF (THICKOPT. EQ. 2) THEN 
DO I=1,NUMPTS 
THICK= THT(I) 

IF (YCORD(I). LT. 0. )YCOR = YCORD( I) -ABS(THICK) 
IF (YCORD(I). GE. 0. )YCOR = YCORD( I)+ABS(THICK) 
YCORD(I) = YCOR 
END DO 
ENDIF 

IF (THICKOPT. EQ. 3) GO TO 64 
GO TO 60 
ENDIF 

CALL CLRSCRN 



OPTION TO MAKE ANOTHER RUN 
PRINT * 

PRINT *, * DO YOU WISH TO: ' 

PRINT *, * 1) MAKE ANOTHER RUN OR* 

PRINT *, * 2) END THIS SESSION* 

PRINT *, * ENTER 1 OR 2. * 

CALL QUERY (NANS) 

IF (NANS .EQ. 1) GO TO 59 
1111 STOP 

920 F0RMAT(1X,F8.4,2E11.4) 

930 FORMAT( 1X,F8. 4,F11. 6) 

1000 FORMAT( /////,* INPUT ALPHA IN DEGREES *) 
1002 FORMAT(Al) 

1006 FORMAT(Il) 

END 






DATA VALUES FOR VARIOUS AIRFOILS. TO USE, REMOVE COMMENTS 
AND PLACE AFTER COMMON CARDS IN MAIN, PROGRAM. 









■ y? Vc it it it it it it it it it it it it it it itititititit 
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